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Abstract. Research efforts of the past fifty years have led to a development of linear 
integer programming as a mature discipline of mathematical optimization. Such a 
level of maturity has not been reached when one considers nonlinear systems subject 
to integrality requirements for the variables. This chapter is dedicated to this topic. 

The primary goal is a study of a simple version of general nonlinear integer 
problems, where all constraints are still linear Our focus is on the computational 
complexity of the problem, which varies significantly with the type of nonlinear 
objective function in combination with the underlying combinatorial structure. Nu- 
merous boundary cases of complexity emerge, which sometimes surprisingly lead 
even to polynomial time algorithms. 

We also cover recent successful approaches for more general classes of problems. 
Though no positive theoretical efficiency results are available, nor are they likely to 
ever be available, these seem to be the currently most successful and interesting 
approaches for solving practical problems. 

It is our belief that the study of algorithms motivated by theoretical considera- 
tions and those motivated by our desire to solve practical instances should and do 
inform one another. So it is with this viewpoint that we present the subject, and it is 
in this direction that we hope to spark further research. 
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1 Overview 

In the past decade, nonlinear integer programming has gained a lot of mindshare. 
Obviously many important applications demand that we be able to handle nonlin- 
ear objective functions and constraints. Traditionally, nonlinear mixed-integer pro- 
grams have been handled in the context of the field of global optimization, where 
the main focus is on numerical algorithms to solve nonlinear continuous optimiza- 
tion problems and where integrality constraints were considered as an afterthought, 
using branch-and-bound over the integer variables. In the past few years, however, 
researchers from the field of integer programming have increasingly studied nonlin- 
ear mixed-integer programs from their point of view. Nevertheless, this is generally 
considered a very young field, and most of the problems and methods are not as 
well-understood or stable as in the case of linear mixed-integer programs. 

Any contemporary review of nonlinear mixed-integer programming will there- 
fore be relatively short-lived. For this reason, our primary focus is on a classification 
of nonlinear mixed-integer problems from the point of view of computational com- 
plexity, presenting theory and algorithms for the efficiently solvable cases. The hope 
is that at least this part of the chapter will still be valuable in the next few decades. 
However, we also cover recent successful approaches for more general classes of 
problems. Though no positive theoretical efficiency results are available — nor are 
they likely to ever be available, these seem to be the currently most successful and 
interesting approaches for solving practical problems. It is our belief that the study 
of algorithms motivated by theoretical considerations and those motivated by our 
desire to solve practical instances should and do inform one another. So it is with 
this viewpoint that we present the subject, and it is in this direction that we hope to 
spark further research. 

Let us however also remark that the selection of the material that we dis- 
cuss in this chapter is subjective. There are topics that some researchers associate 
with "nonlinear integer programming" that are not covered here. Among them are 
pseudo-Boolean optimization, max-cut and quadratic assignment as well as general 
0/1 polynomial programming. There is no doubt that these topics are interesting, 
but, in order to keep this chapter focused, we refrain from going into these topics. 
Instead we refer the interested reader to the references f55l on max-cut, |32| for 
recent advances in general 0/1 polynomial programming, and the excellent surveys 
ll29l on pseudo-Boolean optimization and II103I [341 on the quadratic assignment 
problem. 

A general model of mixed-integer programming could be written as 



max/min 



s.t. 



1 -^n) 

■ ,x„) < 



(1) 



g,n{xi,...,X„) < 
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where , . . . ,g,„ : M" ^ M are arbitrary nonlinear functions. However, in parts 
of the chapter, we study a rather restricted model of nonlinear integer program- 
ming, where the nonlinearity is confined to the objective function, i.e., the following 
model: 

max/min f{xi ,...,x„) 

subject to Ax < b (2) 

xgM"i xZ"2, 

where A is a rational matrix and b is a rational vector. It is clear that this model is 
still NP-hard, and that it is much more expressive and much harder to solve than 
integer linear programs. 

We start out with a few fundamental hardness results that help to get a picture of 
the complexity situation of the problem. 

Even in the pure continuous case, nonlinear optimization is known to be hard. 

Theorem 1. Pure continuous polynomial optimization over polytopes (n2 — 0) in 
varying dimension is NP-hard. Moreover, there does not exist a fully polynomial 
time approximation scheme (FPTAS) (unless P — NPj. 

Indeed the max-cut problem can be modeled as minimizing a quadratic form over 
the cube [—1, 1]", and thus inapproximability follows from a result by Hastad ||65]| . 
On the other hand, pure continuous polynomial optimization problems over poly- 
topes (n2 = 0) can be solved in polynomial time when the dimension ni is fixed. 
This follows from a much more general result on the computational complexity of 
approximating the solutions to general algebraic formulae over the real numbers by 
Renegar 1 1 11 1 . 

However, as soon as we add just two integer variables, we get a hard problem 
again: 

Theorem 2. The problem of minimizing a degree-4 polynomial over the lattice 
points of a convex polygon is NP-hard. 

This is based on the NP-completeness of the problem whether there exists a positive 
integer x <c with x^ = a (mod b); see ll53ll4T1l . 

We also get hardness results that are much worse than just NP-hardness. The neg- 
ative solution of Hubert's tenth problem by Matiyasevich 195. 96 j , based on earlier 
work by Davis, Putnam, and Robinson, implies that nonlinear integer programming 
is incomputable, i.e., there cannot exist any general algorithm. (It is clear that for 
cases where finite bounds for all variables are known, an algorithm trivially ex- 
ists.) Due to a later strengthening of the negative result by Matiyasevich (published 
in |76|), there also cannot exist any such general algorithm for even a small fixed 
number of integer variables; see ||4T1|. 

Theorem 3. The problem of minimizing a linear form over polynomial constraints 
in at most 10 integer variables is not computable by a recursive function. 

Another consequence, as shown by Jeroslow ll75l . is that even integer quadratic 
programming is incomputable. 
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Theorem 4. The problem of minimizing a linear form over quadratic constraints in 
integer variables is not computable by a recursive function. 

How do we get positive complexity results and practical methods? One way is 
to consider subclasses of the objective functions and constraints. First of all, for the 
problem of concave minimization or convex maximization which we study in Sec- 
tion |2] we can make use of the property that optimal solutions can always be found 
on the boundary (actually on the set of vertices) of the convex hull of the feasible 
solutions. On the other hand, as in the pure continuous case, convex minimization, 
which we address in Section [3]l, is much easier to handle, from both theoretical and 
practical viewpoints, than the general case. Next, in Section]?] we study the general 
case of polynomial optimization, as well as practical approaches for handling the 
important case of quadratic functions. Finally, in Section|5] we briefly describe the 
practical approach of global optimization. 

For each of these subclasses covered in sections |2j|5] we discuss positive com- 
plexity results, such as polynomiality results in fixed dimension, if available (Sec- 



tions pT^pT 4.1 1, including some boundary cases of complexity in Sections 2.2 
|3.2| and |5.2| and discuss practical algorithms (Sections |2.3||3.3[|4.2||4.3| [5. 

We end the chapter with conclusions (Section [6|, including a table that summa- 
rizes the complexity situation of the problem (Table[T]). 



2 Convex integer maximization 



2.1 Fixed dimension 



Maximizing a convex function over the integer points in a polytope in fixed dimen- 
sion can be done in polynomial time. To see this, note that the optimal value is 
taken on at a vertex of the convex hull of all feasible integer points. But when the 
dimension is fixed, there is only a polynomial number of vertices, as Cook et al. Il38l 
showed. 

Theorem S. Let P = {x eR" : Ax <h} be a rational polyhedron with A £ Q™^" 
and let (j) be the largest binary encoding size of any of the rows of the system Ax < b. 
Let = conv(Pn Z") be the integer hull of P. Then the number of vertices ofP^ is 
at most2m"{6n^(j))" 

Moreover, Hartmann ll64l gave an algorithm for enumerating all the vertices, which 
runs in polynomial time in fixed dimension. 

By using Hartmann's algorithm, we can therefore compute all the vertices of the 
integer hull P', evaluate the convex objective function on each of them and pick the 
best. This simple method already provides a polynomial-time algorithm. 
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2.2 Boundary cases of complexity 

In the past fifteen years algebraic geometry and commutative algebra tools have 
shown their exciting potential to study problems in integer optimization (see 1191 
11311 and references therein). The presentation in this section is based on the papers 

The first key lemma, extending results of 1 1011 for combinatorial optimization, 
shows that when a suitable geometric condition holds, it is possible to efficiently 
reduce the convex integer maximization problem to the solution of polynomially 
many linear integer programming counterparts. As we will see, this condition holds 
naturally for a broad class of problems in variable dimension. To state this result, we 
need the following terminology. A direction of an edge (i.e., a one-dimensional face) 
e of a polyhedron P is any nonzero scalar multiple of u — v with u, v any two distinct 
points in e. A set of vectors covers all edge-directions ofP if it contains a direction 
of each edge of P. A linear integer programming oracle for matrix A G Z'"^" and 
vector b e Z'" is one that, queried on w € Z", solves the linear integer program 
max{w^x : Ax = b, x g N"}, that is, either returns an optimal solution x e N", or 
asserts that the program is infeasible, or asserts that the objective function w^x is 
unbounded. 

Lemma 1. For any fixed d there is a strongly polynomial oracle-time algorithm 
that, given any vectors Wi, . . . ,Wj G Z", matrix A e J/"^" and vector b G Z™ en- 
dowed with a linear integer programming oracle, finite set E <zTL" covering all 
edge-directions of tlie polyhedron conv{x G N" : Ax = b}, and convex functional 
c : — > M presented by a comparison oracle, solves the convex integer program 

max{c(w|x, . . . ,wjx) : Ax = b, x G N"} . 

Here, solving the program means that the algorithm either returns an optimal 
solution X G N", or asserts the problem is infeasible, or asserts the polyhedron 
{x G Wl : Ax = b} is unbounded; and strongly polynomial oracle-time means 
that the number of arithmetic operations and calls to the oracles are polynomially 
bounded in m and «, and the size of the numbers occurring throughout the algorithm 
is polynomially bounded in the size of the input (which is the number of bits in the 
binary representation of the entries of wi, . . . jW^/jAjbjis). 

Let us outline the main ideas behind the proof to Lemma [T] and let us point 
out the difficulties that one has to overcome. Given data for a convex integer 
maximization problem max{c(wi'^x, . . . , wjx) : Ax = b, x G N"}, consider the 
polyhedron P := conv{x G N" : Ax = b} C M" and its linear transformation Q 
{(wj^x, . . . , wjx) : X G P} C M''. Note that P has typically exponentially many ver- 
tices and is not accessible computationally. Note also that, because c is convex, there 
is an optimal solution x whose image (wj'^x, . . . , wJx) is a vertex of Q. So an impor- 
tant ingredient in the solution is to construct the vertices of Q. Unfortunately, Q may 
also have exponentially many vertices even though it lives in a space M'' of fixed di- 
mension. However, it can be shown that, when the number of edge-directions of P 
is polynomial, the number of vertices of Q is polynomial. Nonetheless, even in this 
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case, it is not possible to construct these vertices directly, because the number of 
vertices of P may still be exponential. This difficulty can finally be overcome by 
using a suitable zonotope. See BSlllOll for more details. 

Let us now apply Lemma [T| to a broad (in fact, universal) class of convex inte- 
ger maximization problems. Lemma[T]implies that these problems can be solved in 
polynomial time. Given an (r + s) x t matrix A, let Ai be its r x f sub-matrix consist- 
ing of the first r rows and let A2 be its s x f sub-matrix consisting of the last s rows. 
Define the n-fold matrix of A to be the following {r + ns) x nt matrix. 



(l««)Ai)®(/„®A2) 



/Ai Ai Ai ••• Ai\ 
A2 ••• 
A2 ••• 

V • • • A2 y 



Note that A'"' depends on r and s: these will be indicated by referring to A as an 
"(r + i) X t matrix." 

We establish the following theorem, which asserts that convex integer maximiza- 
tion over «-fold systems of a fixed matrix A, in variable dimension nt, are solvable 
in polynomial time. 

Theorem 6. For any fixed positive integer d and fixed {r + s) x t integer matrix A 
there is a polynomial oracle-time algorithm that, given n, vectors wi, . . . , G Z"' 

and b G and convex function c : M'' > M presented by a comparison oracle, 

solves the convex n-fold integer maximization problem 

max{c(w|x, . . . , wjx) : A^"'^x = b, x e N'"} . 

The equations defined by an n-fold matrix have the following, perhaps more illu- 
minating, interpretation: splitting the variable vector and the right-hand side vec- 
tor into components of suitable sizes, x = (x^,....x") and b — (b'^jb', . . . ,b"), 
where b" G 1/ and x*^ € N' and b*^ G 1/ for k = !,....«, the equations become 
^\(Sk=i x*^) = b'^ and A2X*^ — b*^ for k — Thus, each component x*^ satis- 

fies a system of constraints defined by A2 with its own right-hand side b*^, and the 
sum L^^jX*^ obeys constraints determined by A\ and b'' restricting the "common 
resources shared by all components". 

Theorem|6]has various applications, including multiway transportation problems, 
packing problems, vector partitioning and clustering. For example, we have the fol- 
lowing corollary providing the first polynomial time solution of convex 3-way trans- 
portation. 

Corollary 1. (convex 3-way transportation) For any fixed d,p,q there is a poly- 
nomial oracle-time algorithm that, given n, arrays wi , . . . , S Z''^'?^", u G Z^^^, 

V € Z''^", z G Tfl^", and convex c : M presented by comparison oracle, solves 

the convex integer 3-way transportation problem 
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max{c(w|x,...,wjx) : xeN'""'''" , ^JC; j,*: = Zi,<r , J^Xij^k^Vi^k , = m,j } ■ 

i j k 

Note that in contrast, if the dimensions of two sides of the tables are variable, 
say, q and n, then already the standard linear integer 3-way transportation problem 
over such tables is NP-hard, see [44, 45 , 46 1. 

In order to prove Theorem [6] we need to recall some definitions. The Graver 
basis of an integer matrix A G Z™^", introduced in ll59l . is a canonical finite set ^#(A) 
that can be defined as follows. For each of the 2" orthants of M" let Hj denote the 
inclusion-minimal Hilbert basis of the pointed rational polyhedral cone ker(A) n 
ffj. Then the Graver basis ^(A) is defined to be the union ^(A) = U^" \ {0} 
over all these 2" Hilbert bases. By this definition, the Graver basis ^^(A) has a nice 
representation property: every z € ker(A) n Z" can be written as a sign-compatible 
nonnegative integer linear combination z = a,g, of Graver basis elements g, e 
^^(A). This follows from the simple observation that z has to belong to some orthant 
i^j of K" and thus it can be represented as a sign-compatible nonnegative integer 
linear combination of elements in Hj C ^#(A). For more details on Graver bases and 
the currently fastest procedure for computing them see 1 124 67l ll30l| . 

Graver bases have another nice property: They contain all edge directions in the 
integer hulls within the polytopes f\, = {x : Ax = b, x > 0} as b is varying. We 
include a direct proof here. 

Lemma 2. For every integer matrix A G U"^" and every integer vector b G N'", the 
Graver basis ^(A) of A covers all edge-directions of the polyhedron convjx e N" : 
Ax = b } defined by A and b. 

Proof. Consider any edge e of P := convjx G N" : Ax = b} and pick two distinct 
points u, V G e n N". Then g := u — v is in ker(A) n Z" \ {0} and hence, by the 
representation property of the Graver basis ?#(A), g can be written as a finite sign- 
compatible sum g — ^g' with g' G ^^(A) for all /. Now, we claim that u — g' G P 
for all /. To see this, note first that g' G J#(A) C ker(A) implies Ag' = and hence 
A(u — g') = Au = b; and second, note that u — g' > 0: indeed, if g'- < then uj — 
g'j > Uj > 0; and if g'j > then sign-compatibility implies g'j < gj and therefore 
Uj-g'j>uj~gj^vj>0. 

Now let w G M" be a linear functional uniquely maximized over P at the edge e. 
Then for all /, as just proved, u — g' G P and hence w^g' > 0. But £ w^g' = w^g = 
w^u — w^v — 0, implying that in fact, for all /, we have w^g' = and therefore 
u — g' G e. This implies that each g' is a direction of the edge e (in fact, moreover, 
all g' are the same, so g is a multiple of some Graver basis element). 



In Section |3.2[ we show that for fixed matrix A, the size of the Graver basis of 
a'"' increases only polynomially in n implying Theorem 14 that states that certain 



integer convex n-fold minimization problems can be solved in polynomial time when 
the matrix A is kept fixed. As a special case, this implies that the integer linear n- 
fold maximization problem can be solved in polynomial time when the matrix A is 
kept fixed. Finally, combining these results with Lemmas[T]and|2j we can now prove 
Theorem |6] 
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Proof (of Theorem^. The algorithm underlying Proposition 14 provides a poly- 
nomial time realization of a linear integer programming oracle for A*^"' and b. The 
algorithm underlying Proposition [2] allows to cornpute the Graver basis ^^(A^")) in 
time which is polynomial in the input. By Lemmaol this set £ ^^(A'"') covers all 
edge-directions of the polyhedron conv{x £ N"' : AJ"'x = b } underlying the convex 
integer program. Thus, the hypothesis of Lemma [T] is satisfied and hence the algo- 
rithm underlying Lemma [T] can be used to solve the convex integer maximization 
problem in polynomial time. 



2.3 Reduction to linear integer programming 

In this section it is our goal to develop a basic understanding about when a discrete 
polynomial programming problem can be tackled with classical linear integer op- 
timization techniques. We begin to study polyhedra related to polynomial integer 
programming. The presented approach applies to problems of the kind 

max{/(x) : Ax = b, X e Z'|} 

with convex polynomial function /, as well as to models such as 

max{c"^x -.xeKi), (3) 

where K[ denotes the set of all integer points of a compact basic-closed semi- 
algebraic set K described by a finite number of polynomial inequalities, i. e., 

K^{x(^W : pi{x) <Q, ieM, 1 < x < u}. 

We assume that /;>,■ S Z[x] :— Z[xi, . . . ,x„], for all / e M = {1, . . . ,m}, and l,u G Z". 

One natural idea is to derive a linear description of the convex hull of Kj. Un- 
fortunately, com{Kj) might contain interior integer points that are not elements of 
K, see Figure [T| On the other hand, if a description of conv(^r/) is at hand, then the 
mathematical optimization problem of solving Q can be reduced to optimizing the 
linear objective function c^x over the polyhedron conv(A7). This is our first topic 
of interest. In what follows we denote for a set D the projection of D to a set of 
variables x by the symbol Dx- 

Definition 1. For polynomials pi , . . . ,/7„, : Z" ^ Z we define the polyhedron asso- 
ciated with the vector p(x) = {pi{x),... ,p,n{x)) of polynomials as 

Pp = conv({(x,p(x)) eM"+'" |xe [i,u]nz"}). 

The significance of the results below is that they allow us to reformulate the non- 
linear integer optimization problem max{/(x) : Ax = b, x e Z']_} as the optimization 
problem minjTT ; Ax = b, /(x) < TT, x e Z'|}. This in turn has the same objective 
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Fig. 1 Let a: = {x G M- I xiX2 - 1 < 0, < ati < 3, < X2 < 5}. Tlie point (1,2) is contained in 
conv(A7). But it violates the constraint xiX2 — 1 < 0. 



function value as the linear integer program: min{7r :Ax~b, (x, n) E Pf, x G 

In this situation, an H-description or V-description of the polyhedron Pf is sufficient 

to reduce the original nonlinear optimization problem to a linear integer problem. 

Proposition 1. For a vector of polynomials p e Z[x]"', let 



Proof. It suffices to show that Kj C (Pp n {(x, It) £ M"+'" : It < 0})^. Let us consider 
X e X/ C [l,u] nZ". By definition, (x,p(x)) e Pp. Moreover, we have p(x) < 0. This 
impHes (x, p(x)) e (Pp n {(x, it) G W+'" : TT < 0}) , and thus, 

xe (Ppn{(x,;r)eM"+"':;r<0})^. 

Note that even in the univariate case, equality in Formula Q of Proposition [T] 
does not always hold. For instance, if Kj — {x £ Z : — 5 < 0, — 3 < x < 5}, then 

com{Ki) = [-2,2] ^ [-2.2,2.2] C (p,.2 n {(x, ;/r) G : - 5 < 0})^. 

Although even in very simple cases the sets com{Kj) and (Pp n {(x, ;r) : ;r < 0})^ 
differ, it is still possible that the integer points in Kj and (Pp fl {(x, It) : Jt < 0})^ are 
equal. In our example with = {x G Z : — 5 < 0, — 3 < x < 5}, we then obtain. 



Kj = {xeZ" : p(x) < 0, 1 < X < u}. 



Then, 




(4) 



iT/ = {-2,-1,0, 1,2} = [-2.2,2.2] nZ. 



Of course, for any p G Z[x]'" we have that 



10 Raymond Hemmecke, Matthias Koppe, Jon Lee and Robert Weismantel 

/:/ = {xeZ" :p(x) <o, i<x<u} c (Ppn{(x,;r) : ;r<0})^nz". (5) 

The key question here is when equahty holds in Formula (j5]). 

Theorem 7. Let p e Z[x]'" and Kj^{xe Z" : p(x) < 0, 1 < x < u}. Then, 

Ki = (Ppn{(x,;r):;r<0})^ n Z" 

holds if every polynomial p' e {pi : i = 1, . . . ,m} satisfies the following condition 
y(EAkk) - EAky(k) < 1, (6) 

k k 

forallXk>0, ke [l,u]nZ«, IkAk = 1 and Ik^kkGZ". 

Proof Using Formula (jsj, we have to show that (fp n {(x, Tt) : Tt <0})^nZ" C K[ 
if all Pi, i e {!,..., m}, satisfy Equation Let x e (fpn{(x,;r) G M"+" : ;r < 
0})j^nZ". Then, there exists alt e M.'" such that 

(x, ;r) e Pp n {(x, it) e : ;r < 0}. 

By definition, ;r < 0. Furthermore, there must exist nonnegative real numbers Ak > 
0, k e [l,u] nZ", such that Y^kh. = 1 and (x,lt) = £k'^k(k,p(k)). Suppose that 
there exists an index io such that the inequality pi^ (x) < is violated. The fact that 
Pig € Z[x] and x G Z", implies that Pigix) > I. Thus, we obtain 

J^hpioi^) = TTio < < 1 < p,g{x) = Pig (^Akk) , 

k k 

or equivalently, Pig {'Yk^kk) —Hk^kPioik) > 1. Because this is a contradiction to 
our assumption, we have that p,(x) < for all /. Hence, x G Kj. This proves the 
claim. 

The next example illustrates the statement of Theorem|7j 

Example 1. Let p G Z[x], x i— > p(x) 2>x\ +2x^ — 19. We consider the semi- 
algebraic set 

{xG I /7(x) < 0, <xi < 3, <ji:2 < 3} and Ki = KnZ^. 

It turns out that the convex hull of Kj is described by xi +X2 < 3, < xi < 2 and 
Q<X2- Notice that the polynomial p is convex. This condition ensures that p satisfies 
Equation (|6| of Theorem|7] We obtain in this case 

{Ppn{{x,7:) em.^ : n <0})^^ {xgM^ ; 9x1 +6x2 < 29, 

3x1 + 10x2 < 31, 
9x1 + 10x2 < 37, 
15x1 +2x2 < 37, 
15x1+6x2 <41, 
-XI < 0, <X2 < 3}. 
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The sets K, com{K[) and [Pp n {(x, tt) e : tt < 0})^ ai'e illustrated in Figure[2] 
Note that here = {P,, C^{{^,n) : n <Q})^C^ 1?. 




Next we introduce a large class of nonhnear functions for which one can ensure 
that equality holds in Formula (jSj. 

Definition 2. {Integer-Convex Polynomials) 

A polynomial / e Z[x] is said to be integer-convex on [l,u] C M", if for any finite 
subset of nonnegative real numbers {Ak}kG[i,u]nZ" ^ ^+ with LkG[i,u]nZ" ■^k = 1 and 
Ek6[i,u]nz" ■^kk e [l,u] nZ", the following inequality holds: 

/( E Akk) < ^k/(k). (7) 

kG[i,u]nZ'" kG[i,u]nZ'" 

If (|7ji holds strictly for all {Ak}kG[i,u]nZ" ^ andx e [l,u] nZ" such that ^k^k — 
1, X = £k Akk, and Ax < 1, then the polynomial / is called strictly integer-convex on 
[l,u]. 

By definition, a (strictly) convex polynomial is (strictly) integer-convex. Con- 
versely, a (strictly) integer-convex polynomial is not necessarily (strictly) convex. 
Figure[3]gives an example. 

Integer convexity is inherited under taking conic combinations and applying a 
composition rule. 

(a) For any finite number of integer-convex polynomials e Z [x] , s £ { 1 , . . . , f } , on 
[1, u] , and nonnegative integers Os G Z+ , s G {1, . . . ,t}, the polynomial / G Z [x] , 
X 1-^ /(x) :— Y!s=i is integer-convex on [l,u]. 
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Fig. 3 The graph of an integer-convex polynomial p on [1,4]. 



(b) Letl,u e Z" andlet e Z[x], x i-> /j(x) := c^x + y, be a linear function. Setting 
"W = {c^x + 7 : x G [1, u] }, for every integer-convex univariate polynomial q G 
Z[w], the function p E Z[x], x p{x) := q{h{x)) is integer-convex on [l,u]. 

Indeed, integer-convex polynomial functions capture a lot of combinatorial struc- 
ture. In particular, we can characterize the set of all vertices in an associated poly- 
hedron. Most importantly, if / is integer-convex on [l,u], then this ensures that for 
any integer point x e [I, u] the value /(x) is not underestimated by all TT G M with 
(x,7r) G Pf, where Pf is the polytope associated with the graph of the polynomial 

/ez[x]. 

Theorem 8. For a polynomial f G Z[x] and l,u G Z", 1 -|- 1 < u, let 

Pf = conv ({ (x,/(x)) G Z"+i I X G [l,u] nZ«}) . 

Then, f is integer-convex on [1, u] is equivalent to the condition that for all (x, n) G 
Pf, X G Z" we have that f{x) < n. Moreover, if f is strictly integer-convex on [l,u], 
then for every x G [l,u] flZ", the point (x,/(x)) is a vertex ofPf. 

Proof. First let us assume that / is integer-convex on [1, u] . Let (x, 7t) G Pf such that 
X G Z". Then, there exist nonnegative real numbers {Ak}k6[i,u]nz" ^ ILk^k = 1, 
such that (x, tt) = Y.k ^k (k, /(k) ) . It follows that 

/W = /(lAkk) < £Ak/(k) = 7t. 

^ k ^ k 

Next we assume that / is not integer-convex on [I,u]. Then, there exists a subset 
of normegative real numbers {Ak}kg[i_u]nz" ^ ^+ witii Dk^k = 1 such that 
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X := ^Akk e [l,u] nZ" and % := ^Ak/(k) < /(^Akk) = /(x). 

k k ^ k ^ 

But then, (x,7r) = Lk'^k(k,/(k)) G P/ violates the inequality /(x) < %. This is a 
contradiction to the assumption. 

If / is strictly integer-convex on [1, u], then for each x G [1, u] n Z", we have that 

/(x)< Y. ^k/(k), 

kG[l,u]nZ"\{x} 

for all Ak e M+, k e [l,u] nZ" \ {x}, with x = Ik^kk and Lk'^-k = 1- Thus, every 
point (x,/(x)), X e [l,u] nZ", is a vertex of Pf. 

3 Convex integer minimization 

The complexity of the case of convex integer minimization is set apart from the gen- 
eral case of integer polynomial optimization by the existence of bounding results for 
the coordinates of optimal solutions. Once a finite bound can be computed, it is clear 
that an algorithm for minimization exists. Thus the fundamental incomputability re- 
sult for integer polynomial optimization (Theorem[3]l does not apply to the case of 
convex integer minimization. 

The first bounds for the optimal integer solutions to convex minimization prob- 
lems were proved by fTS! '1251 ■ We present the sharpened bound that was obtained 
by ifTTlfTOl for the more general case of quasi-convex polynomials. This bound is a 
consequence of an efficient theory of quantifier elimination over the reals; see 111 lOL 

Theorem 9. Let f,gi,...,g,„€Z [xi,. . . ,Xn] be quasi-convex polynomials of degree 
at most d>2, whose coefficients have a binary encoding length of at most £. Let 

F= {xeM" :g,(x) <0 fori^\,...,m) 

be the (continuous) feasible region. If the integer minimization problem min{/(x) : 
X g F n Z" } is bounded, there exists a radius R G Z+ of binary encoding length at 
most {md)'^^"^l such that 

min{/(x) :xeFnZ"} =min{/(x) :xeFnZ", ||x|| <R]. 



3.1 Fixed dimension 

In fixed dimension, the problem of convex integer minimization can be solved us- 
ing variants of Lenstra's algorithm ifSTl for integer programming. Indeed, when the 
dimension n is fixed, the bound R given by Theorem |9] has a binary encoding size 
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that is bounded polynomially by the input data. Thus, a Lenstra-type algorithm can 
be started with a "small" (polynomial-size) initial outer ellipsoid that includes a 
bounded part of the feasible region containing an optimal integer solution. 

The first algorithm of this kind for convex integer minimization was announced 
by Khachiyan |78|. In the following we present the variant of Lenstra's algorithm 
due to Heinz |66|, which seems to yield the best known complexity bound for the 
problem. The complexity result is the following. 

Theorem 10. Let f,g\,... ,^,„ G l^xy ,.. . ^x„] be quasi-convex polynomials of de- 
gree at most d >2, whose coefficients have a binary encoding length of at most L 
There exists an algorithm running in time m£^^^^ d^^"'^2'^^"^ ^ that computes a mini- 
mizer x* e Z" of the problem ([TJ or reports that no minimizer exists. If the algorithm 
outputs a minimizer li*, its binary encoding size is £d^'^"\ 

A complexity result of greater generality was presented by Khachiyan and Porko- 
lab fW|. It covers the case of minimization of convex polynomials over the inte- 
ger points in convex semialgebraic sets given by arbitrary (not necessarily quasi- 
convex) polynomials. 

Theorem 11. Let Y (-R^ be a convex set given by 

F = { y e : Qjx' e M"' : • • • Q^x" e R"'^ : P(y,x' , . . . ,x'*') } 

with quantifiers Qj G {3,V}, where P is a Boolean combination of polynomial in- 
equalities 

g,(y,xS...,x'«) <0, 

with degrees at most d > 2 and coefficients of binary encoding size at most L 
There exists an algorithm for solving the problem minjy^^ : y G FflZ*} in time 

£0(i)M)0('^')nr=,o(«,). 

When the dimension k-\-YJiL\ni is fixed, the algorithm runs in polynomial time. 
For the case of convex minimization where the feasible region is described by 



convex polynomials, the complexity bound of Theorem 1 1 however, translates to 



£0(i),„o(«2)^o(«'')^ which is worse than the bound of Theorem 10 |66 



In the remainder of this subsection, we describe the ingredients of the variant 
of Lenstra's algorithm due to Heinz. The algorithm starts out by "rounding" the 
feasible region, by applying the shallow-cut ellipsoid method to find proportional 
inscribed and circumscribed ellipsoids. It is well-known |60| that the shallow-cut el- 
lipsoid method only needs an initial circumscribed ellipsoid that is "small enough" 
(of polynomial binary encoding size - this follows from Theorem |9]l and an imple- 
mentation of a shallow separation oracle, which we describe below. 

For a positive-definite matrix A we denote by <^(A,x) the ellipsoid {x e M" : 
(x-x)^A(x-x) < 1}. 

Lemma 3 (Shallow separation oracle). Let gQ,...,g,„+\ £ Z[x] be quasi-convex 
polynomials of degree at most d, the binary encoding sizes of whose coefficients are 
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Fig. 4 The implementation of the shallow separation oracle, (a) Test points x, y in the circumscribed 
ball cf (1,0). (b) Case I: All test points x,] are (continuously) feasible; so their convex hull (a cross- 
polytope) and its inscribed ball S((n+ 1)^^^,0) are contained in the (continuous) feasible region F . 



at most r. Let the ( continuous) feasible region F = {x eM." : gi(x) <0} be contained 
in the ellipsoid (o^A^x), where A and x have binary encoding size at most L There 
exists an algorithm with running time m(lnr)'^^^^ d'^^"^ that outputs 

(a) "true " if 

£'{{n+l)-^A,x)CF C^{A,x); (8) 

(b) otherwise, a vector c G Q" \ {0} of binary encoding length (Z + r){dn)'^^^^ with 

f C(?(A,x)n{xeM":c"^(x-x)< ^(c"^Ac)^/2|_ (9) 

Proof. We give a simplified sketch of the proof, without hard complexity estimates. 
By applying an affine transformation to F C S'{A,x), we can assume that F is con- 
tained in the unit ball J'(/,0). Let us denote as usual by ei,...,e„ the unit vec- 
tors and by e„+i,...,e2„ their negatives. The algorithm first constructs numbers 
1,1, • ■ -Aid > with 

3 < An < • • • < A,y < — — (10) 

n+§ n+1 

and the corresponding point sets B, = {x,y :— XjjCi : j ~ 1, . . . ,d}; see Figure|4](a). 
The choice of the bounds ( fTO] ) for X^j will ensure that we either find a large enough 
inscribed ball for (a) or a deep enough cut for (b). Then the algorithm determines 
the (continuous) feasibility of the center and the 2n innermost points x, j . 

Case I. If X, 1 e F for / = l,...,2n, then the cross-polytope conv{x, i : / — 
1, . . . ,2n} is contained in F; see Figure |4](b). An easy calculation shows that the 
ball (?((«+ 1)^^,0) is contained in the cross-polytope and thus in F; see Figure [4] 
Hence the condition in (a) is satisfied and the algorithm outputs "true". 

Case II. We now discuss the case when the center violates a polynomial in- 
equality gdix) < (say). Let = {x G M" : go{x) <0}^F. Due to convexity of Fq, 
for all / = 1 , . . . , n, one set of each pair B, n Fo and B„+, n Fq must be empty; see Fig- 
ure |5]( a). Without loss of generality, let us assume Bn+jOFo = for all /. We can 
determine whether a «-variate polynomial function of known maximum degree d is 
constant by evaluating it on (d+l)" suitable points (this is a consequence of the 
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Fig. 5 The implementation of the shallow separation oracle, (a) Case II: The center violates a 
polynomial inequality go (x) < (say). Due to convexity, for all i = 1 , . . . , n, one set of each pair 
Bi r\F and B„+i HF must be empty, (b) Case III: A test point xi^[ is infeasible, as it violates an 
inequality go(x) < (say). However, the center is feasible at least for this inequality. 



Fundamental Theorem of Algebra). For our case of quasi-convex polynomials, this 
can be improved; indeed, it suffices to test whether the gradient Vgo vanishes on the 
nd points in the set Bi U • • • UB„. If it does, we know that go is constant, thus F — 
and so can we return an arbitrary vector c. Otherwise, there is a point x,j e B, with 
c :— Vgo(xjy) 7^ 0; we return this vector as the desired normal vector of a shallow 
cut. Due to the choice of A,y as a number smaller than the cut is deep enough 
into the ellipsoid S{A,x), so that (|9]l holds. 

Case III. The remaining case to discuss is when G F but there exists a e 
{1, . . . ,2n} withx^ 1 ^ F. Without loss of generality, let ^= l,andletxi_i violate the 
polynomial inequality go(x) < 0, i.e., go(xi,i ) > 0; see Figurel5](b). We consider the 
univariate polynomial <^{X) = gQ{Xti). We have 0(0) = gQ{0)< and 0(Ai.i) > 0, 
so ([) is not constant. Because (/) has degree at most d, its derivative 0' has degree at 
most d — I, so has at most d — 1 roots. Thus, for at least one of the d different 
values . . ,Ai^d, say we must have 0'(Aij) 7^ 0. This implies that c 
V^o(xi y) 0. By convexity, we have xij ^ F, so we can use c as the normal vector 
of a shallow cut. 

By using this oracle in the shallow-cut ellipsoid method, one obtains the follow- 
ing result. 

Corollary 2. Let go, . . . ,gm G ^ [x] be quasi-convex polynomials of degree at most d > 
2. Let the ( continuous) feasible region F = { x G M" : g,(x) <0} be contained in the 
ellipsoid ff(Ao,0), given by the positive-definite matrix Aq G Q"^". Let £ G Q>o be 
given. Let the entries of Aq and the coefficients of all monomials of go, . . . ,g,„ have 
binary encoding size at most i. 

There exists an algorithm with running time m{in)^^^^d^^"^ that computes a 
positive-definite matrix A G Q"^" and a point x G Q" with 

(a) either S'{{n+ l)-^A,x) C F C S'{A,x) 

(b) orF C <f (A,x) andvolS'{A,x) < e. 

Finally, there is a lower bound for the volume of a continuous feasible region F 
that can contain an integer point. 
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Lemma 4. Under the assumptions of Theorem^ if F H Z" ^ 0, then there exists an 
e £ Q>o of binary encoding size i{dn)°^^^ with volF > e. 

On the basis of these results, one obtains a Lenstra-type algorithm for the de- 
cision version of the convex integer minimization problem with the desired com- 
plexity. By applying binary search, the optimization problem can be solved, which 



provides a proof of Theorem 10 



3.2 Boundary cases of complexity 

In this section we present an optimality certificate for problems of the form 

min{/(x) :Ax = b,l<x<u,xeZ"}, 
where A e Z''^", b G Z'', l,u G Z", and where / : M" ^ M is a separable convex 

n 

function, that is, /(x) = Y, fii^i) with convex functions /; : E ^ M, / = 1, . . . ,«. 

1=1 

This certificate then immediately leads us to a oracle-polynomial time algorithm to 
solve the separable convex integer minimization problem at hand. Applied to sep- 
arable convex n-fold integer minimization problems, this gives a polynomial time 
algorithm for their solution |69|. 

For the construction of the optimality certificate, we exploit a nice super-additivity 
property of separable convex functions. 

Lemma 5. Let f : W ^ R be a separable convex function and let hi ,.. .hi, £ W 
belong to a common orthant ofW, that is, they all have the same sign pattern from 
{> 0, < 0}". Then, for any x G M" w have 




Proof. The claim is easy to show for n = 1 by induction. If, w.l.o.g., h2 > h\ > 
then convexity of /implies [f{x + h\+h2) — f{x + h2)]lh\ > [f{x + h\) — f{x)\/h\, 
and thus /(jc + Zzi +/12) -/(x) > [/(x + /i2) -/(x)] + [f{x + hx)-f{x)\. The claim 
for general n then follows from the separability of / by adding the superadditivity 
relations of each one-parametric convex summand of /. 

A crucial role in the following theorem is again played by the Graver basis J^(A) 
of A. Let us remind the reader that the Graver basis J^(A) has a nice representa- 
tion property due to its definition: every z G ker(A) nZ" can be written as a sign- 
compatible nonnegative integer linear combination z — ^, a,g, of Graver basis el- 
ements g; G '^{A). This followed from the simple observation that z has to belong 
to some orthant ffj of W and thus it can be represented as a sign-compatible non- 
negative integer linear combination of elements in Hj C J#(A) belonging to this 
orthant. Note that by the integer Caratheodory property of Hilbert bases, at most 
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2 • dim(ker(A)) — 2 vectors are needed in such a representation II118L It is precisely 
this simple representation property of $^(A) combined with the superadditivity of 
the separable convex function / that turns '^(A) into an optimality certificate for 
min{/(x) : Ax = b,l < X < u,x e Z"}. 

Theorem 12. Let / : M" ^ M /je a separable convex function given by a comparison 
oracle that when queried on x,y Cz Z" decides whether /(x) < /(y), /(x) = /(y), 
or f(x) > /(y). Then xq is an optimal feasible solution to min{/(x) : Ax = b,l < 
X < u,x G Z"} if and only if for all g € ^^(A) the vector Xq + g is not feasible or 
/(xo + g)>/(xo). 

Proof. Assume that xq is not optimal and let x^in be an optimal solution to the given 
problem. Then Xmjn — xq G ker(A) and thus it can be written as a sign-compatible 
nonnegative integer linear combination Xmin — xq = (^iSi of Graver basis elements 
gi S ^^(A). We show that one of the g, must be an improving vector, that is, for some 
gi we have that xq + g/ is feasible and /(xq + g,) < /(xq). 

For all /, the vector g, has the same sign-pattern as x^in — xq and it is now easy to 
check that the coordinates of xq + gi lie between the corresponding coordinates of 
xo and Xmm- This implies in particular 1 < xq + g, < u. Because g, £ ker(A), we also 
have A(xo + g,) = b for all Consequently, for all / the vector xq + g, would be a 
feasible solution. It remains to show that one of these vectors has a strictly smaller 
objective value than xq. 

Due to the superadditivity from Lemma |5j we have 

C2n-2 \ k 
X0+ E «'&■) -/(xo)>i;«'[/(xo + g/)-/(xo)]. 

Thus, at least one of the summands /(xq + g,) — /(xq) must be negative and we have 
found an improving vector for zq in ^#(A). 

We now turn this optimality certificate into a polynomial oracle-time algorithm 
to solve the separable convex integer minimization problem min{/(x) : Ax = b,l < 
X < u,x G Z"}. For this, we call ag a greedy Graver improving vector if xq + ag is 
feasible and such that /(xq + ag) is minimal among all such choices of a G Z+ and 
g G ?#(A). Then the following result holds. 

Theorem 13. Let f :M." ^Rbe a separable convex function given by a comparison 
oracle. Moreover, assume that |/(x)| < M for all x £ {x : Ax = b,l < x < u,x G 
Z"}. Then any feasible solution xq to min{/(x) : Ax = b,l < x < u,x G Z"} can be 
augmented to optimality by a number of greedy Graver augmentation steps that is 
polynomially bounded in the encoding lengths of A, b, 1, u, M, and xq. 

Proof. Assume that xo is not optimal and let Xmin be an optimal solution to the given 
problem. Then Xmin — xq G ker(A) and thus it can be written as a sign-compatible 
nonnegative integer linear combination x^in — xq = ^, a,g, of at most 2n — 2 Graver 



basis elements g,- G ^^(A). As in the proof of Theorem 12 sign-compatibility implies 



that for all / the coordinates of xq + o:,g, lie between the corresponding coordinates 
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of xo and Xmin- Consequently, we have 1 < xq + a,g, < u. Because g,- G ker(A), we 
also have A(xo + 0£;g,) = b for all /. Consequently, for all / the vector xq + 0£,g,- would 
be a feasible solution. 

Due to the superadditivity from Lemma |5] we have 

(2«-2 \ k 
X0+ L a,g, 1 -/(xo) > ^^[/(xo + a,g,)-/(xo)]. 

Thus, at least one of the summands /(xq + o:,g,) — /(xq) must be smaller than 
2;^[/(xmin) " /(xo)], giving an improvement that is at least 2^Z2 times the max- 
imal possible improvement /(xmin) — /(xq). Such a geometric improvement, how- 
ever, implies that the optimum is reached in a number of greedy augmentation steps 
which is polynomial in the encoding lengths of A, b, 1, u, M, and xq 13. 

Thus, once we have a polynomial size Graver basis, we get a polynomial time 
algorithm to solve the convex integer minimization problem at hand. 

For this, let us consider again n-fold systems (introduced in Section [Z2| . Two 
nice stabilization results established by Ho§ten and Sullivant [711 and Santos and 
Sturmfels 11131 immediately imply that if Ai and A2 are kept fixed, then the size of 
the Graver basis increases only polynomially in the number n of copies of A 1 and 
A2. 

Proposition 2. For any fixed {r + s) x t integer matrix A there is a polynomial time 
algorithm that, given any n, computes the Graver basis 'i^(A^"^) of the n-fold matrix 
AW = (l„0Ai)e(/„(g)A2). 



Combining this proposition with Theorem 13 we get following nice result from 



Theorem 14. Let A be a fixed integer (r + s) x t matrix and let f : R"' R be 
any separable convex function given by a comparison oracle. Then there is a poly- 
nomial time algorithm that, given n, a right-hand side vector b G Z''+"' and some 
bound |/(x) \ < M on f over the feasible region, solves the n-fold convex integer 
programming problem 

min{/(x) : A^x^^b, xeN'"}. 

Note that by applying an approach similar to Phase 1 of the simplex method one 
can also compute an initial feasible solution xq to the n-fold integer program in 
polynomial time based on greedy Graver basis directions P3ll67l . 

We wish to point out that the presented approach can be generalized to the mixed- 
integer situation and also to more general objective functions that satisfy a certain 
superadditivity/subadditivity condition, see ll68l[86l for more details. Note that for 
mixed-integer convex problems one may only expect an approximation result, as 
there need not exist a rational optimum. In fact, already a mixed-integer greedy aug- 
mentation vector can be computed only approximately. Nonetheless, the technical 
difficulties when adjusting the proofs for the pure integer case to the mixed-integer 
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situation can be overcome ll68l . It should be noted, however, that the Graver basis 
of n-fold matrices does not show a stabihty result similar to the pure integer case 
as presented in ITTTIII 1311 . Thus, we do not get a nice polynomial time algorithm for 
solving mixed-integer convex n-fold problems. 



3.3 Practical algorithms 



In this section, the methods that we look at, aimed at formulations having convex 
continuous relaxations, are driven by O.R./engineering approaches, transporting and 
motivated by successful mixed-integer linear programming technology and smooth 



continuous nonlinear programming technology. In Section 3.3.1 we discuss general 
algorithms that make few assumptions beyond those that are typical for convex con- 



tinuous nonlinear programming. In Section 3.3.2 we present some more specialized 
techniques aimed at convex quadratics. 



3.3.1 General algorithms 

Practical, broadly applicable approaches to general mixed-integer nonlinear pro- 
grams are aimed at problems involving convex minimization over a convex set with 
some additional integrality restriction. Additionally, for the sake of obtaining well- 
behaved continuous relaxations, a certain amount of smoothness is usually assumed. 
Thus, in this section, the model that we focus on is 

min /(x,y) 

f^'l:' (p[i."i) 

i<y <u 

xeM"', yeZ"2, 

where / : M" ^ M and g: M" ^ M™ are twice continuously-differentiable convex 
functions, 1 e (ZU {-°o})"2, u G (ZU {+°o})«2^ and 1 < u. It is also helpful to 
assume that the feasible region of the relaxation of (P[l,u]) obtained by replacing 
y e Z"- with y e M"2 is bounded. We denote this continuous relaxation by (Ph[1,u]). 

To describe the various algorithmic approaches, it is helpful to define some re- 
lated subproblems of (P[l,u]) and associated relaxations. Our notation is already 
designed for this. For vector 1' e (ZU{-°o})"2 and u' G (ZU{+°o})"2, with 
1 < 1' < u' < u, we have the subproblem (P[l',u']) and its associated continuous 
relaxation (Pr[1',u']). 

Akeady, we can see how the family of relaxations (Pr[1',u']) leads to the ob- 
vious extension of the Branch-and-Bound Algorithm of mixed-integer linear pro- 
gramming. Indeed, this approach was experimented with in [63 1. The Branch-and- 
Bound Algorithm for mixed-integer nonlinear programming has been implemented 
as MINLP-BB ll88ll . with continuous nonlinear-programming subproblem relax- 
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ations solved with the active-set solver f ilterSQP and also as SBB, with asso- 
ciated subproblems solved with any of CONOPT, SNOPT and MINOS. Moreover, 
Branch-and-Bound is available as an algorithmic option in the actively developed 
code Bonmin 121112511231 . which can be used as a callable library, as a stand-alone 
solver, via the modeling languages AMPL and GAMS, and is available in source-code 
form, under the Common Public License, from COIN-OR |28 1, available for running 
on NEOS |27 1. By default, relaxations of subproblems are solved with the interior- 
point solver Ipopt (whose availability options include all of those for Bonmin), 
though there is also an interface to filters QP. The Branch-and-Bound Algorithm 
in Bonmin includes effective strong branching and SOS branching. It can also be 
used as a robust heuristic on problems for which the relaxation (Pr) does not have 
a convex feasible region, by setting negative 'cutoff gaps'. 

Another type of algorithmic approach emphasizes continuous nonlinear pro- 
gramming over just the continuous variables of the formulation. For fixed y e Z"^, 
we define 

min /(x,y) 
s.t. g(x,y) < 

y = y 
xeM"i. 

Clearly any feasible solution to such a continuous nonlinear-programming subprob- 
lem (P'') yields an upper bound on the optimal value of (P[l,u]). When (P^) is 
infeasible, we may consider the continuous nonlinear-programming feasibility sub- 
problem 

m 

min E 

1=1 

s.t. g(x,y)<w - 

y = y 

If we can find a way to couple the solution of upper-bounding problems (P^) (and 
the closely related feasibility subproblems (F^)) with a lower-bounding procedure 
exploiting the convexity assumptions, then we can hope to build an iterative proce- 
dure that will converge to a global optimum of (P[l, u]). Indeed, such a procedure is 
the Outer-Approximation (OA) Algorithm Il48ll49l . Toward this end, for a finite set 
of "linearization points" 

|(x^ eM"i,y*^eM"2^ : ;t = 1, . . . , 



we define the mixed-integer linear programming relaxation 
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mm z 

s.t. v/(x^/')T {^-''^^ +/(x^y*) < z, v (x^/■) e 

Vg(x^y^)T +g(x^y^) < 0, V (x^y*) G (P'^[l,u]) 

X e M"i 

yeM"2, l<y<u 

zeM. 

We are now able to concisely state the basic OA Algorithm. 
Algorithm 1 (OA Algorithm) 

Input: The mixed-integer nonlinear progratn (P[l,u]). 
Output: An optimal solution (x*,y*). 

1. Solve the nonlinear-programming relaxation (-Pr), let (x',y') be an optimal 
solution, and let K := 1, so that initially we have = {{x\y^)}. 

2. Solve the mixed-integer linear programming relaxation (f'^[l,u]), andlet {n* ,y* ,z*) 
be an optimal solution. If (x*,y*,z*) corresponds to a feasible solution of 
(P[l,u]) (i.e, iff(x*,y*) < z* and g(x*,y*) < Oj, then STOP (with the optimal 
solution (x*,y*) o/(P[l,u])j. 

3. Solve the continuous nonlinear-programming subproblem [P^ ). 

i. Either a feasible solution (x*,y*,z*) is obtained, 

ii. or [P^ ) is infeasible, in which case we solve the nonlinear-programming 
feasibility subproblem (F^), and let its solution be (x*,y*,u*) 

4. In either case, we augment the set of linearization points, by letting K := 
K+land (x^,y^) :=(x*,y*). 

5. GOTO^ 



Each iteration of Steps [3]|4] generate a linear cut that can improve the mixed- 
integer linear programming relaxation (P*^ [1, u] ) that is repeatedly solved in Step 
[2] So clearly the sequence of optimal objective values for (P'^[l,u]) obtained in 
Step |2] corresponds to a nondecreasing sequence of lower bounds on the optimum 
value of (P[l,u]). Moreover each linear cut returned from Steps 3]|4 cuts off the 



previous solution of (P''^ [1, u] ) from Step [2] A precise proof of convergence (see 
for example |21|) uses these simple observations, but it also requires an additional 
assumption that is standard in continuous nonlinear programming (i.e. a "constraint 
qualification"). 

Implementations of OA include D I COPT ll39l which can be used with either of 
the mixed-integer linear programs codes Cplex and Xpress-MP, in conjunction 
with any of the continuous nonlinear programming codes CONOPT, SNOPT and 
MINOS and is available with GAMS. Additionally Bonmin has OA as an algorith- 
mic option, which can use Cplex or the COIN-OR code Cbc as its mixed-integer 
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linear programming solver, and Ipopt or FilterSQP as its continuous nonlinear 
programming solver. 

Generalized Benders Decomposition f54\ is a technique that is closely related 
to and substantially predates the OA Algorithm. In fact, one can regard OA as a 
proper strengthening of Generalized Benders Decomposition (see Il48ll49l ). so as a 
practical tool, we view it as superseded by OA. 

Substantially postdating the development of the OA Algorithm is the simpler 
and closely related Extended Cutting Plane (ECP) Algorithm introduced in II135L 
The original ECP Algorithm is a straightforward generalization of Kelley's Cutting- 
Plane Algorithm ||77ll for convex continuous nonlinear programming (which pre- 
dates the development of the OA Algorithm). Subsequently, the ECP Algorithm has 
been enhanced and further developed (see, for example II 1361 11341 ) to handle, for 
example, even pseudo-convex functions. 

The motivation for the ECP Algorithm is that continuous nonlinear programs are 
expensive to solve, and all that the associated solutions give us are further lineariza- 
tion points for (P'^[l,u]). So the ECP Algorithm dispenses altogether with the so- 
lution of continuous nonlinear programs. Rather, in the most rudimentary version, 
after each solution of the mixed-integer linear program (P'^[l,u]), the most vio- 
lated constraint (i.e, of /(x*,y*) < z and §(x*,y*) < 0) is linearized and appended 
to (P^[l,u]). This simple iteration is enough to easily establish convergence (see 
111351 ). It should be noted that for the case in which there are no integer-constrained 
variables, then at each step (P^[l,u]) is just a continuous linear program and we 
exactly recover Kelley's Cutting-Plane Algorithm for convex continuous nonlinear 
programming. 

It is interesting to note that Kelley, in his seminal paper llTTll . already considered 
application of his approach to integer nonlinear programs. In fact, Kelley cited Go- 
mory's seminal work on integer programming |'58^,'57 | which was also taking place 
in the same time period, and he discussed how the approaches could be integrated. 

Of course, many practical improvements can be made to the rudimentary ECP 
Algorithm. For example, more constraints can be linearized at each iteration. An im- 
plementation of the ECP Algorithm is the code Alpha-ECP (see 113411 ) which uses 
Cplex as its mixed-integer linear programming solver and is available with GAMS. 
The general experience is that for mildly nonlinear problems, an ECP Algorithm can 
outperform an OA Algorithm. But for a highly nonlinear problem, the performance 
of the ECP Algorithm is limited by the performance of Kelley's Cutting-Plane Al- 
gorithm, which can be quite poor on highly-nonlinear purely continuous problems. 
In such cases, it is typically better to use an OA Algorithm, which will handle the 
nonlinearity in a more sophisticated manner. 

In considering again the performance of an OA Algorithm on a mixed-integer 
nonlinear program (P[l,u]), rather than the convex continuous nonlinear program- 
ming problems (P^ ) and (F'*) being too time consuming to solve (which led us 
to the ECP Algorithm), it can be the case that solution of the mixed-integer linear 
programming problems (P''*^[l,u]) dominate the running time. Such a situation led 
to the Quesada-Grossmann Branch-and-Cut Algorithm I1107t . The viewpoint is that 
the mixed-integer linear programming problems iV'^ [1, u] ) are solved by a Branch- 
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and-Bound or Branch-and-Cut Algorithm. During the solution of the mixed-integer 
linear programming problem (P^ [1, u] ), whenever a new solution is found (i.e., one 
that has the variables y integer), we interrupt the solution process for (P^[l,u]), and 
we solve the convex continuous nonlinear programming problems (P^ ) to derive 
new outer approximation cuts that are appended to mixed-integer linear program- 
ming problem (P'^ [I7 u] ) . We then continue with the solution process for (P^ [1, u] ) . 
The Quesada-Grossmann Branch-and-Cut Algorithm is available as an option in 
Bonmin. 

Finally, it is clear that the essential scheme of the Quesada-Grossmann Branch- 
and-Cut Algorithm admits enormous flexibility. The Hybrid Algorithm \2l \ incor- 
porates two important enhancements. 

First, we can seek to further improve the linearization (P^[l,u]) by solving con- 
vex continuous nonlinear programming problems at additional nodes of the mixed- 
integer linear programming Branch-and-Cut tree for (P ''^ [1, u] ) — that is, not just 
when solutions are found having y integer In particular, at any node (P'^[l',u']) of 
the mixed-integer linear programming Branch-and-Cut tree, we can solve the as- 
sociated convex continuous nonlinear programming subproblem (Pk[1',u']): Then, 
if in the solution (x* , y* ) we have that y* is integer, we may update the incumbent 
and fathom the node; otherwise, we append (x*,y*) to the set of linearization 
points. In the extreme case, if we solve these continuous nonlinear programming 
subproblems at every node, we essentially have the Branch-and-Bound Algorithm 
for mixed-integer nonlinear programming. 

A second enhancement is based on working harder to find a solution (x* , y* ) with 
y* integer at selected nodes (P^[l',u']) of the mixed-integer linear programming 
Branch-and-Cut tree. The idea is that at a node (P^[l',u']), we perform a time- 
limited mixed-integer linear programming Branch-and-Bound Algorithm. If we are 
successful, then we will have found a solution to the node with (x*.y*) with y* 
integer, and then we perform an OA iteration (i.e.. Steps 3]|4i on (P[l',u']) which 
will improve the linearization (P''^[l',u']). We can then repeat this until we have 
solved the mixed-integer nonlinear program (P[l',u']) associated with the node. If 
we do this without time limit at the root node (P[l,u]), then the entire procedure 
reduces to the OA Algorithm. The Hybrid Algorithm was developed for and first 
made available as part of Bonmin. 

FilMint fD is another successful modern code, also based on enhancing the 
general framework of the Quesada-Grossmann Branch-and-Cut Algorithm. The 
main additional innovation introduced with FilMint is the idea of using ECP 
cuts rather than only performing OA iterations for getting cuts to improve the 
linearizations (P'^[l',u']). Subsequently, this feature was also added to Bonmin. 
FilMint was put together from the continuous nonlinear programming active-set 
code FilterSQP, and the mixed-integer linear programming code MINTO. Cur- 
rently, FilMint is only generally available via NEOS |[50l . 

It is worth mentioning that just as for mixed-integer linear programming, effec- 
tive heuristics can and should be used to provide good upper bounds quickly. This 
can markedly improve the performance of any of the algorithms described above. 
Some examples of work in this direction are |24 | and [221 . 
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3.3.2 Convex quadratics and second-order cone programming 

Though we will not go into any details, there is considerable algorithmic work and 
associated software that seeks to leverage more specialized (but still rather gen- 
eral and powerful) nonlinear models and existing convex continuous nonlinear- 
programming algorithms for the associated relaxations. In this direction, recent 
work has focused on conic programming relaxations (in particular, the semi -definite 
and second-order cones). On the software side, we point to work on the binary 
quadratic and max-cut problems (via semi-definite relaxation) 11081 1109 1 with the 
code Biq Mac EOl . We also note that Cplex (vll) has a capability aimed at 
solving mixed-integer quadratically-constrained programs that have a convex con- 
tinuous relaxation. 

One important direction for approaching quadratic models is at the modeling 
level. This is particulary useful for the convex case, where there is a strong and ap- 
pealing relationship between quadratically constrained programming and second- 
order cone programming (SOCP). A second-order cone constraint is one that ex- 
presses that the Euclidean norm of an affine function should be no more than another 
affine function. An SOCP problem consists of minimizing a linear function over a 
finite set of second-order cone constraints. Our interest in SOCP stems from the fact 
that (continuous) convex quadratically constrained programming problems can be 
reformulated as SOCP problems (see |92|). The appeal is that very efficient interior- 
point algorithms have been developed for solving SOCP problems (see ||56l , for 
example), and there is considerable mature software available that has functional- 
ity for efficient handling of SOCP problems; see, for example: SDPT3 ||117J (GNU 
GPL open-source license; Matlab) , SeDuMi II119I (GNU GPL open-source license; 
Matlab), LOQO |93 | (proprietary; C library with interfaces to AMPL and Matlab), 
MOSEK L99| (proprietary; C library with interface to Matlab), Cplex Ell (propri- 
etary; C library). Note also that MOSEK and Cplex can handle integer variables as 
well; one can expect that the approaches essentially marry specialized SOCP solvers 
with Branch-and-Bound and/or Outer-Approximation Algorithms. Further branch- 
and-cut methods for mixed-integer SOCP, employing linear and convex quadratic 
cuts 1 36 1 and a careful treatment of the non-differentiability inherent in the SOCP 
constraints, have recently been proposed pTl]. 

Also in this vein is recent work by Giinliik and Linderoth ||6T1 |62l . Among 
other things, they demonstrated that many practical mixed-integer quadratically 
constrained programming formulations have substructures that admit extended for- 
mulations that can be easily strengthened as certain integer SOCP problems. This 
approach is well known in the mixed-integer linear programming literature. Let 



n 



vv e M, X e M';, z e {0, 1}" : w > £ 




> Xi > liZi, / = 1,2, . . . ,n 



!=1 



where r, g M+ and m,-,/, g M for all / = 1,2, ... ,n. The set Q appears in several 
formulations as a substructure. Consider the following extended formulation of Q 
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{xi,yi,Zi) £ Si, 1,2,...,«|, 



where 



Si := ^{xi,yi,Zi) e X {0, 1} : > xj, mzi > Xi > Uzi, Xi > o|, 
and M, , /, £ M. The convex hull of each 5, has the form 

Si := {{xi,yi,Zi) £ : > xf, UiZi > Xi > Uzi, 1 > Z/ > 0, > O} 

(see ||35ll6Tll62lll29l ). Note that x? —jiZi is not a convex function, but nonetheless 
SI is a convex set. Finally, we can state the result of 1611 l62l . which also follows 
from a more general result of [70], that the convex hull of the extended formulation 
Q has the form 



ye: 



Z £ . 



w > 



ixi,yi,Zi) £ Si, ! = 1,2,...,«|. 



Note that all of the nonlinear constraints describing the SI and Q'^ are rather simple 
quadratic constraints. Generally, it is well known that even the "restricted hyperbolic 
constraint" 

n 

yiZi>Y,xl,x£M",yi>0, zi>0 

(more general than the nonconvexity in SI) can be reformulated as the second-order 
cone constraint 

2x 



yi - Zi 



<yi+Zi 



In this subsection, in the interest of concreteness and brevity, we have focused 
our attention on convex quadratics and second-order cone programming. However, 
it should be noted that a related approach, with broader applicability (to all con- 
vex objective functions) is presented in fSTJ, and a computational comparison is 
available in |52|. Also, it is relevant that many convex non-quadratic functions are 
representable as second-order cone programs (see Bl fTSll ). 



4 Polynomial optimization 

In this section, we focus our attention on the study of optimization models involving 
polynomials only, but without any assumptions on convexity or concavity. It is worth 
emphasizing the fundamental result of Jeroslow (Theorem|4]i that even pure integer 
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quadratically constrained programming is undecidable . One can however avoid this 
daunting pitfall by bounding the variables, and this is in fact an assumption that we 
should and do make for the purpose of designing practical approaches. From the 
point of view of most applications that we are aware of, this is a very reasonable 
assumption. We must be cognizant of the fact that the geometry of even quadratics 
on boxes is daunting from the point of view of mathematical programming; see 
Figure [6] Specifically, the convex envelope of the graph of the product x\X2 on a 
box deviates badly from the graph, so relaxation-based methods are intrinsically 
handicapped. It is easy to see, for example, that for 5i , ^2 > 0, x\X2 is strictly convex 
on the line segment joining (0,0) and ^2); while xiX2 is strictly concave on the 
line segment joining (5i,0) and (0,52) ■ 




Fig. 6 Tetrahedral convex envelope of the graph of the product x\X2 on a box 



Despite these difficulties, we have positive results. In Section 4.1 the highlight 
is ci fully polynomial time approximation scheme (FPTAS) for problems involving 
maximization of a polynomial in fixed dimension, over the mixed-integer points in 
a polytope. In Section |4.2| we broaden our focus to allow feasible regions defined 
by inequalities in polynomials (i.e., semi-algebraic sets). In this setting, we do not 
present (nor could we expect) complexity results as strong as for linear constraints, 
but rather we show how tools of semi-definite programming are being developed to 



provide, in a systematic manner, strengthened relaxations. Finally, in Section 4.3 



we describe recent computational advances in the special case of semi-algebraic 
programming for which all of the functions are quadratic — i.e., mixed-integer 
quadratically constrained programming (MIQCP). 



4.1 Fixed dimension and linear constraints: An FPTAS 



As we pointed out in the introduction (Theorem |2]l, optimizing degree-4 polynomi- 
als over problems with two integer variables is akeady a hard problem. Thus, even 
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when we fix the dimension, we cannot get a polynomial-time algorithm for solving 
the optimization problem. The best we can hope for, even when the number of both 
the continuous and the integer variables is fixed, is an approximation result. 

Definition 3. (FPTAS) 

(a) An algorithm is an approximation algorithm for a maximization problem 
with optimal cost /max, if for each instance of the problem of encoding length n, 

runs in polynomial time in n and returns a feasible solution with cost f^, 

such that /c/ > ( 1 - e) • /max- 

(b) A family {£4}e of £-approximation algorithms is 2. fully polynomial time ap- 
proximation scheme (FPTAS) if the running time of si/^ is polynomial in the 
encoding size of the instance and l/e. 

Indeed it is possible to obtain an FPTAS for general polynomial optimization of 
mixed-integer feasible sets in polytopes 141114011421 . To explain the method of the 
FPTAS, we need to review the theory of short rational generating functions pio- 
neered by Barvinok lfl2HT3l . The FPTAS itself appears in Section ]?. 1.3| 

4.1.1 Introduction to rational generating functions 

We explain the theory on a simple, one-dimensional example. Let us consider the 
set S of integers in the interval f = [0, . . . , n] ; see the top of Figure|7](a). We associate 

with S the polynomial g{S;z) — + ~\ hz"^' i.e., every integer a £ S 

corresponds to a monomial z" with coefficient 1 in the polynomial g{S;z)- This 
polynomial is called the generating function of S (or of P). From the viewpoint 
of computational complexity, this generating function is of exponential size (in the 
encoding length of n), just as an explicit list of all the integers 0, l,...,n— l,n 
would be. However, we can observe that g{S;z) is a finite geometric series, so there 
exists a simple summation formula that expresses it in a much more compact way: 

g(5;z) = z" + z' + --- + z"~'+z"= ^7^"^' . (11) 

1 — z 

The "long" polynomial has a "short" representation as a rational function. The en- 
coding length of this new formula is linear in the encoding length of n. 

Suppose now someone presents to us a finite set S of integers as a generating 
function g{S;z)- Can we decide whether the set is nonempty? In fact, we can do 
something much stronger even - we can count the integers in the set S, simply by 

evaluating at g{S;z) at z = 1. On our example we have |5|=g(5;l) = l'' + l'H h 

1"^' + 1" = n + 1. We can do the same on the shorter, rational-function formula if 
we are careful with the (removable) singularity z = 1. We just compute the limit 

l-z"+' -fn+l)z" 
151 = \img(S;z) = lim — = lim — ^ =n+l 

z^l z^l 1-Z z^l -1 
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Fig. 7 (a) One-dimensional Brion theorem, (b) The domains of convergence of the Laurent series. 



using the Bernoulli-FHopital rule. Note that we have avoided to carry out a poly- 
nomial division, which would have given us the long polynomial again. 

The summation formula ( [TT| i can also be written in a slightly different way: 

1 z"+' 1 z" 

8iS;z) = - ^ = ^ + T^ (12) 

1— z 1-z 1— z 1— z' 

Each of the two summands on the right-hand side can be viewed as the summation 
formula of an infinite geometric series: 

gi{z) = j^=z'' + z'+z^ + ..., (13a) 
i — z 



82{z) = = z" +z"-' +z"-2 . . . . (13b) 

i — z 



The two summands have a geometrical interpretation. If we view each geometric 
series as the generating function of an (infinite) lattice point set, we arrive at the 
picture shown in Figure [T] We remark that all integer points in the interval [0,n] 
are covered twice, and also all integer points outside the interval are covered once. 
This phenomenon is due to the one-to-many correspondence of rational functions 
to their Laurent series. When we consider Laurent series of the function (z) about 
z = 0, the pole z = 1 splits the complex plane into two domains of convergence 
(Figure [rji: For |z| < 1, the power series z" + z' +z^ + . . . converges to giiz). As a 
matter of fact, it converges absolutely and uniformly on every compact subset of the 
open circle {z S C : |z| < 1 }. For |z| > 1, however, the series diverges. On the other 
hand, the Laurent series — z^' —z^^ — z^^ — . . . converges (absolutely and compact- 
uniformly) on the open circular ring {z G C : |z| > 1 } to the function g\ (z), whereas 
it diverges for |z| < 1. The same holds for g2{z)- Altogether we have: 



(14) 
(15) 
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(a) 




(b) 
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Fig. 8 (a) Tiling a rational two-dimensional cone with copies of the fundamental parallelepiped, 
(b) The semigroup SCI? generated by bi and b2 is a linear image of 



We can now see that the phenomenon we observed in formula ( [T3] l and Figure |7] 
is due to the fact that we had picked two Laurent series for the summands gi{z) 
and g2{z) that do not have a common domain of convergence; the situation of for- 
mula ^T3\ and Figure|7]appears again in the cZ-dimensional case as Brion 's Theorem. 

Let us now consider a two-dimensional cone C spanned by the vectors bi ~ 
(a, — 1) and b2 — 1); see Figure [s] for an example with a — 2 and j3 = 4. We 
would like to write down a generating function for the integer points in this cone. We 
apparently need a generalization of the geometric series, of which we made use in 
the one-dimensional case. The key observation now is that using copies of the half- 
open fundamental parallelepiped, n = { Aibi + A2b2 : Ai e [0, l),X2 G [0, 1) }, the 
cone can be tiled: 

C=\J{s + n) where 5 = { ^iibi + Ai2b2 : (ah ,Ai2) G } (16) 

seS 

(a disjoint union). Because we have chosen integral generators bi,b2, the integer 
points are "the same" in each copy of the fundamental parallelepiped. Therefore, 
also the integer points of C can be tiled by copies of 17 n Z^; on the other hand, we 
can see Cnl? as a finite disjoint union of copies of S, shifted by the integer points 
ofH; 

cnz^ = y (s+(nnz2)) = |J (x+s). (17) 

The set S is just the image of under the matrix (bi,b2) = (^j ^); cf. Figure 
[s] Now is the direct product of Z+ with itself, whose generating function is 
the geometric series g(Z+;z) = + + + • • • = j^. We thus obtain the 

generating function as a product, g{Z^;zi,Z2) = j^j^ • Applying the linear 
transformation (bi,b2). 
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(a) (b) "> (c) "1 

Fig. 9 (a) A cone of index 5 generated by b' and b^. (b) A triangulation of the cone into the 
two cones spanned by {b',w} and {b^,w}, having an index of 2 and 3, respectively. We have 
the inclusion-exclusion formula g(cone{bi ,b2};z) = g(cone{bi,w};z) + g(cone{b2,w};z) — 
g(cone{w};z); here the one-dimensional cone spanned by w needed to be subtracted, (c) A 
signed decomposition into the two unimodular cones spanned by {b' , w'} and {b^, w'}. We have 
the inclusion-exclusion formula g(cone{bi,b2};z) = g(cone{bi,w'};z) — g(cone{b2, w'};z) + 
g(cone{w'};z). 



From ([17} it is now clear that ^(C;zi,Z2) = Hxennz^ ^7 8{S;zi ,Z2); the multipli- 
cation with the monomial z^'z,^ corresponds to the shifting of the set S by the vector 
{xi ,X2). In our example, it is easy to see that HnZ^ = { (/,0) : / = 0, . . . , a-f j3 — 1 }. 
Thus 

(l-zfZ2')(l-zfzl) 

Unfortunately, this formula has an exponential size as the numerator contains a + P 
summands. To make the formula shorter, we need to recursively break the cone into 
"smaller" cones, each of which have a much shorter formula. We have observed 
that the length of the formula is determined by the number of integer points in the 
fundamental parallelepiped, the index of the cone. Triangulations usually do not 
help to reduce the index significantly, as another two-dimensional example shows. 
Consider the cone C' generated by b i = ( 1 , 0) and b2 = ( 1 , o: ) ; see Figure |9] We have 
n' r\l? — { (0, 0) } U { ( 1 , /) : i — 1 , . . . , a — 1 } , so the rational generating function 
would have a summands in the numerator, and thus have exponential size. Every 
attempt to use triangulations to reduce the size of the formula fails in this example. 
The choice of an interior vector w in Figure|9j for instance, splits the cone of index 5 
into two cones of index 2 and 3, respectively - and also a one-dimensional cone. 
Indeed, every possible triangulation of C' into unimodular cones contains at least a 
two-dimensional cones! The important new idea by Barvinok was to use so-called 
signed decompositions in addition to triangulations in order to reduce the index of 
a cone. In our example, we can choose the vector w = (0, 1) from the outside of the 
cone to define cones Ci = cone{bi,w} and C2 — cone{w,b2}; see Figure[9] Using 
these cones, we have the inclusion-exclusion formula 
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(a) 



(b) 




Fig. 10 Brion's theorem, expressing the generating function of a polyhedron to those of the sup- 
porting cones of all vertices 

g{C';zi,Z2) ^ g{Ci;zi,Z2) - g{C2;zi,Z2) + g{CinC2;zi,Z2) 

It turns out that all cones Ci and C2 are unimodular, and we obtain the rational 
generating function by summing up those of the subcones, 



4.1.2 Barvinok's algorithm for short rational generating functions 

We now present the general definitions and results. Let P C be a rational 
polyhedron. We first define its generating function as the formal Laurent series 
g{P;z) =LaG/>nZ''^" ^ '^[[zi,---,Zd,Zi^,...,Zj^]], i.e., without any consideration 
of convergence properties. (A formal power series is not enough because monomi- 
als with negative exponents can appear) As we remarked above, this encoding of a 
set of lattice points does not give an immediate benefit in terms of complexity. We 
will get short formulas only when we can identify the Laurent series with certain 
rational functions. Now if P is a polytope, then giP^z) is a Laurent polynomial (i.e., 
a finite sum of monomials with positive or negative integer exponents), so it can be 
naturally identified with a rational function g{P;z). Convergence comes into play 
whenever P is not bounded, since then g{P\z) can be an infinite formal sum. We 
first consider a pointed polyhedron P, i.e., P does not contain a straight line. 

Theorem 15. Let P C be a pointed rational polyhedron. Then there exists a non- 
empty open subset U C C' such that the series g{P; z) converges absolutely and uni- 
formly on every compact subset of U to a rational function g{P; z) G Q(zi , . . . , z^;). 

Finally, when P contains an integer point and also a straight line, there does not 
exist any point z e where the series g{P',i) converges absolutely. In this case 
we set g{P; z) = 0; this turns out to be a consistent choice (making the map P 1-^ 
g{P\z) a valuation, i.e., a finitely additive measure). The rational function g{P;z) G 
Q(zi, . . . defined as described above is called the rational generating function 

off nZ''. 



g{C';zi,Z2) 



1 1 1 



(1-Zl)(l-Z2) (1-Z1Z2«)(1-Z2) l~ZlZr 
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Theorem 16 (Brion ||30I)- l^^t P be a rational polyhedron and V{P) be the set of 
vertices of P. Then, g{P;z) — Lvgv(/')^(C/'(v);z), where C/>(v) = v + cone(P — v) 



is the supporting cone of the vertex v; see Figure 10 



We remark that in the case of a non-pointed polyhedron P, i.e., a polyhedron that 
has no vertices because it contains a straight line, both sides of the equation are zero. 

Barvinok's algorithm computes the rational generating function of a polyhe- 
dron P as follows. By Brion's theorem, the rational generating function of a polyhe- 
dron can be expressed as the sum of the rational generating functions of the support- 
ing cones of its vertices. Every supporting cone v, + C, can be triangulated to obtain 
simplicial cones v, + C,y. If the dimension is fixed, these polyhedral computations 
all run in polynomial time. 

Now let K be one of these simplicial cones, whose basis vectors bi,...,bf/ 
(i.e., primitive representatives of its extreme rays) are the columns of some ma- 
trix B e Z''^''; then the index of K is |detB|. Barvinok's algorithm now computes a 
signed decomposition of K to produce simplicial cones with smaller index. To this 

end, it constructs a vector w = aibi H h Udhd G If' \ {0} with |a, | < |detB|"^/''. 

The existence of such a vector follows from Minkowski's first theorem, and it can 
be constructed in polynomial time using integer programming or lattice basis re- 
duction followed by enumeration. The cone is then decomposed into cones spanned 
by d vectors from the set {bi, . . . ,bd,w}; each of the resulting cones then has an 
index at most (ind/T)'''^ In general, these cones form a signed decomposition 
of K; only if w lies inside K, they form a triangulation (see Figure|9]l. The resulting 
cones and their intersecting proper faces (arising in an inclusion-exclusion formula) 
are recursively processed, until cones of low index (for instance unimodular cones) 
are obtained. Finally, for a unimodular cone v + B]R'{_, the rational generating func- 
tion is z''/nf=i(l ~'^')^ where a is the unique integer point in the fundamental 
parallelepiped. We summarize Barvinok's algorithm below. 

Algorithm 2 (Barvinok's algorithm) 

Input: A polyhedron P C M'' given by rational inequalities. 
Output; The rational generating function for P {^1/ in the form 

where £,■ G {±1}, a,- G Z'', and b// G H' . 

1. Compute all vertices v,- and corresponding supporting cones Ci of P. 

2. Triangulate Ci into simplicial cones Cij, keeping track of all the intersecting 
proper faces. 

3. Apply signed decomposition to the cones V; + C,/ to obtain unimodular cones 
V; + Ciji, keeping track of all the intersecting proper faces. 

4. Compute the unique integer point a, in the fundamental parallelepiped of every 
resulting cone v; + Ciji. 

5. Write down the formula ( |18[ ). 
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Fig. 11 Approximation properties of ^^-norms 



We remark that it is possible to avoid computations with the intersecting proper faces 
of cones (step 2 of the algorithm) entirely, using techniques such as polarization, 
irrational decomposition lISOl . or half-open decomposition lISTIISn . 

Due to the descent of the indices in the signed decomposition procedure, the 
depth of the decomposition tree is at most 1 1 + I where D = IdetBl. Be- 

cause at each decomposition step at most 0(2^) cones are created and the depth of 
the tree is doubly logarithmic in the index of the input cone, Barvinok could obtain 
a polynomiality result infixed dimension: 

Theorem 17 (Barvinok |12|). Let d be fixed. There exists a polynomial-time algo- 
rithm for computing the rational generating function \\%) of a polyhedron P C W' 
given by rational inequalities. 



4,1.3 The FPTAS for polynomial optimization 

We now describe the fully polynomial-time approximation scheme, which appeared 
in Il4ni40ll42ll . It makes use of the elementary relation 



max{ii,...,iA,}= lim ^4' (19) 

which holds for any finite set 5 = {ii , . . . , i^f} of non-negative real numbers. This 
relation can be viewed as an approximation result for ^^^-norms. Now if P is a poly- 
tope and / is an objective function non-negative on f n I/' , let x' , . . . , denote all 
the feasible integer solutions in P n If' and collect their objective function values 
Si = /(x') in a vector s G Q^. Then, comparing the unit balls of the 4-norm and the 



£oo-norm (Figure 1 1 1, we get the relation 

L,,:^N-'l''\\sU <UU <Mk=-Uk. 

Thus, for obtaining a good approximation of the maximum, it suffices to solve a 
summation problem of the polynomial function h — f'' on f nZ'' for a value of k that 
is large enough. Indeed, for A: = [(1 + 1 /e) logA^] , we obtain Uk — La- < e/(x™'"') . 
On the other hand, this choice of k is polynomial in the input size (because 1 /e is 
encoded in unary in the input, and logA^ is bounded by a polynomial in the binary 
encoding size of the polytope P). Hence, when the dimension d is fixed, we can 
expand the polynomial function as a list of monomials in polynomial time. 
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Solving the summation problem can be accomplished using short rational gener- 
ating functions as follows. Let g{P',z) be the rational generating function of PnZ'', 
computed using Barvinok's algorithm. By symbolically applying differential op- 
erators to g{P;z), we can compute a short rational function representation of the 
Laurent polynomial g{P,h;z) ~ Y.aePni.'' ^i^)'^'^ ' where each monomial z" corre- 
sponding to an integer point OC G PdZ'' has a coefficient that is the value . To 
illustrate this, consider again the generating function of the interval P = [0,4], 

l—Z i—z 

We now apply the differential operator z^ and obtain 

Applying the same differential operator again, we obtain 

%) [%)"'^'^ = l^^+4z^ + 9z^ + 16/= 

We have thus evaluated the monomial function h{a) = for a = ,4; the 

results appear as the coefficients of the respective monomials. The same works for 
several variables, using the partial differential operators z; ^ for j = 1 , . . . , J. In fixed 
dimension, the size of the rational function expressions occuring in the symbolic 
calculation can be bounded polynomially. Thus one obtains the following result. 

Theorem 18. (a) Let h{xi ,Xd) = ]L/5 cjjx^ G Q[xi , . . . ,Xd] be a polynomial. De- 
fine the differential operator 

Then D/, maps the generating function g(P;z) = T^aePnZ'' '^"^ weighted 
generating function {Dhg){z) =g{P,h;z) = ^^^^pp^rf /z(a)z". 
(b) Let the dimension d be fixed. Let g{P;z) be the Barvinok representation of the 
generating function T^aePnz'' ^" ^f^ ^ G Q [xi , . . . , x^] fee a polynomial, 

given as a list of monomials with rational coefficients cp encoded in binary and 
exponents j3 encoded in unary. We can compute in polynomial time a Barvinok 
representation g{P,h;z)for the weighted generating function Loemz<'''(''^)^"- 

Thus, we can implement the following algorithm in polynomial time. 

Algorithm 3 (Computation of bounds for the optimal value) 

Input: A rational convex polytope P C M''; a polynomial objective function f €E 
Q[xi,...,x^] that is non-negative over f flZ^, given as a list of monomials with 
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rational coefficients cp encoded in binary and exponents /3 encoded in unary; an 
index k, encoded in unary. 

Output: A lower bound Lj^ and an upper bound U^for the maximal function value 
f* of f over Pn Z"^. The bounds Lj^form a nondecreasing, the bounds Uk a nonin- 
creasing sequence of bounds that both reach f* in a finite number of steps. 



1 . Compute a short rational function expression for the generating function g(P\ z) ~ 
^jjgpp^rfz". Using residue techniques, compute |PnZ''| ^ g{P;l) from g{P;z). 

2. Compute the polynomial f^ from f. 

3. From the rational function g{P',i) compute the rational function representation 
ofg{Pj'-z) of^„^p^^,f'{a)z'' by Theorem 
compute 



18 



Using residue techniques. 



g{PJk-\)/g{P-l) 



and Uk '■= 



Taking the discussion of the convergence of the bounds into consideration, one ob- 
tains the following result. 

Theorem 19 (Fully polynomial-time approximation scheme). Let the dimension 
d be fixed. Let P C M'' be a rational convex polytope. Let f be a polynomial with 
rational coefficients that is non-negative on PC] TL'^ , given as a list of monomials 
with rational coefficients cp encoded in binary and exponents /3 encoded in unary. 

(i) Algorithm^computes the bounds Lk, Uk in time polynomial in k, the input size 
ofP and f, and the total degree D. The bounds satisfy the following inequality: 



Uk-Lk<f* 



ipnZ' 



(ii) For k = {1 + l/e)log(|f nZ^'|) (a number bounded by a polynomial in the 
input size), L^ is a (1 — e) -approximation to the optimal value f* and it can 
be computed in time polynomial in the input size, the total degree D, and 1 /e. 
Similarly, Uk gives a (1 + £)-approximation to f*. 

(iii) With the same complexity, by iterated bisection ofP, we can also find a feasible 
solution Xe'EPHZ'' with |/(x£) - /* | < £/* . 

The mixed-integer case can be handled by discretization of the continuous vari- 
ables. We illustrate on an example that one needs to be careful to pick a sequence of 
discretizations that actually converges. Consider the mixed-integer linear optimiza- 
tion problem depicted in Figure 1 2 whose feasible region consists of the point ( 5 , 1 ) 
and the segment { (x,0) : x e [0, 1] }. The unique optimal solution is x = ^, z = I. 
Now consider the sequence of grid approximations where x e ^Z>o. For even m, 
the unique optimal solution to the grid approximation is x= \,z= 1 . However, for 
odd m, the unique optimal solution is x = 0, z = 0. Thus the full sequence of the 
optimal solutions to the grid approximations does not converge because it has two 



limit points; see Figure 12 
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Fig. 12 A mixed-integer linear optimization problem and a sequence of optimal solutions to grid 
problems with two limit points, for even m and for odd m 

To handle polynomial objective functions that take arbitrary (positive and nega- 
tive) values on the feasible region, one can shift the objective function by a constant 
that is large enough. Then, to obtain a strong approximation result, one iteratively 
reduces the constant by a factor. Altogether we have the following result. 

Theorem 20 (Fully polynomial-time approximation schemes). Let the dimension 
n — ni+n2 be fixed. Let an optimization problem (|2| of a polynomial fiinction f over 
the mixed-integer points of a polytope P and an error bound £ be given, where 

(11) / is given as a list of monomials with rational coefficients c^ encoded in 
binary and exponents J3 encoded in unary, 

(12) P is given by rational inequalities in binary encoding, 

(13) the rational number ^ is given in unary encoding. 

(a) There exists a fully polynomial time approximation scheme (FPTAS) for the max- 
imization problem for all polynomial functions /(x,z) that are non-negative on 
the feasible region. That is, there exists a polynomial-time algorithm that, given 
the above data, computes a feasible solution {xg,Zg) E PC) (M"' x Z"2j with 

|_/(Xf,Z£) /"(Xmaxi Zmax) I — ^fi^max^^max)- 

(b) There exists a polynomial-time algorithm that, given the above data, computes 
a feasible solution (x^jZe) g fn (K"' x Z"^) with 

|/(Xe, Zf) /(Xmax: Zniax) I ^ £ |/'(Xmaxi ^max) /(Xmim ^min) I ■ 



4.2 Semi-algebraic sets and SOS programming 

In this section we use results from algebraic geometry over the reals to provide a 
convergent (and in the case of binary optimization, finite) sequence of semi -definite 
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relaxations for the general polynomial optimization problem over semi-algebraic 
sets. 

Z* = minimize /(x) 

s.t. >0, (•= l,...,m, (20) 

where f,gi e M[x] are polynomials defined as: 

/(x)= ^ /«x«, gi{x)= £ g,ax", 

where there are only finitely many nonzero coefficients /« and gi a- Moreover, let 
K = {x £ M." \ gi(x) > 0, / = 1, . . . ,ot} denote the set of feasible solutions. Note 



that problem ( 20 1 can model binary optimization, by taking /(x) = c x, and taking 
as the polynomials gi{x), ajx — b,, —xj and — +Xj (to model xj —Xj = 0). 
Problem ( [20| can also model bounded integer optimization (using for example the 
equation {xj — lj){xj ~ Ij + I) ■ . . . ■ {xj — uj) = to model < xj < uj), as well as 



bounded mixed-integer nonlinear optimization. Problem ( |20| ) can be written as: 
maximize y 

(21) 

s.t. /(x)-7>0, yxeK. 

This leads us to consider conditions for polynomials to be nonnegative over a set K. 

Definition 4. Let p G M[x] where x = (xi , . . . ,x„)^. The polynomial p is called sos 
(sum of squares), if there exist polynomials /ii ,...,/;<. G ]R[x] such that p = Yd^i hj. 

Clearly, in multiple dimensions if a polynomial can be written as a sum of squares 
of other polynomials, then it is nonnegative. However, is it possible for a polyno- 
mial in higher dimensions to be nonnegative without being a sum of squares? The 
answer is yes. The most well-known example is probably the Motzkin-polynomial 
M{x,y,z) = x^y^ +xy + — Sx^z^, which is nonnegative without being a sum 
of squares of polynomials. 

The following theorem establishes a certificate of positivity of a polynomial on 
the set K, under a certain assumption on K. 

Theorem 21 ( II 1061 ll74l ). Suppose that the set K is compact and there exists a poly- 
nomial h(x) of the form 

m 

h{x) = ho{x) + Y^hi{x)gi{x), 

7=1 

such that {x G M" I /i(x) > 0} is compact and /i, (x), / = 0, 1 , . . . , m, are polynomials 
that have a sum of squares representation. Then, if the polynomial g is strictly posi- 
tive over K, then there exist pi G IIi[x], / = 0, 1 , . . . , m, that are sums of squares such 
that 

m 

^(x)=po(x) + £/?,(x)^,(x). (22) 

7=1 
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Note that the number of terms in Equation p2| is Hnear. While the assumption 



of Theorem 21 may seem restrictive, it is satisfied in several cases: 

(a) For binary optimization problems, that is, when K includes the inequalities xj > 



Xj and Xj > xj for all 7=1 



,n. 



(b) If all the gj's are linear, i.e.. A" is a polyhedron. 

(c) If there is one polynomial gi^ such that the set {x G M." \ gk{x) > 0} is compact. 



More generally, one way to ensure that the assumption of Theorem 21 holds 
is to add to K the extra quadratic constraint g„,+i(x) = — ||x|p > for some a 
sufficiently large. It is also important to emphasize that we do not assume that K is 
convex. Notice that it may even be disconnected. 

Let us now investigate algorithmically when a polynomial is a sum of squares. 
As we will see this question is strongly connected to semi-definite optimization. 
The idea of using semi-definite optimization for solving optimization problems over 
polynomials is due to II 1201 and further expanded in ll82l and 110411 . We consider the 
vector 

Vf/(x) — (x ) |or|<(/ — {l,Xl, . . . ,X„,Xi,XlX2, ■ ■ ■ ,X„_lJC„,JC^, . . . ,X[ , . . . ) , 

of all the monomials x" of degree less than or equal to d, which has dimension 

Proposition 3 ( I.37J ). The polynomial g(x) of degree 2d has a sum of squares de- 
composition if and only if there exists a positive semi-definite matrix Q for which 

g(x)=Vrf(x)Tgv,/(x). 

Proof. Suppose there exists an x 5 matrix 2^0 for which ^(x) ~ Vf/(x)^2vrf(x). 
Then Q = HH^ for some s x k matrix H, and thus. 



g(x) = y,ixfHH^y,ix) = £(//Tv,(x))2, 



k 

I 

(=1 



Because {H V£/(x)), is a polynomial, then g(x) is expressed as a sum of squares of 
the polynomials {H~^\ii{x))i. 

Conversely, suppose that g{x) has a sum of squares decomposition g{x) = 
Yfi=ihi{x)'^ . Let h, be the vector of coefficients of the polynomial /i,(x), i.e., 
/z,(x)=hTv,/(x).Thus, 

e 

gix) = £ V,/(x)Th,-hTv./(x) = Vrf(x)Tevj(x), 

1=1 

with Q — jfi^i hihj y 0, and the proposition follows. □ 

Proposition[3]gives rise to an algorithm. Given a polynomial /(x) € K[xi , . . . ,x„] 
of degree 2d. In order to compute the minimum value /* = min{/(x) | x g M"} we 
introduce an artificial variable X and determine 
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max{A I A eM, /(x)-A > 0}. 

With the developments above, we realize that we can determine a lower bound for 
/* by computing the value 

p"" = max{A I A e M, /(x) - A is sos} < /*. 

The latter task can be accomplished by setting up a semi-definite program. In fact, 
if we denote by /« the coefficient of the monomial x" in the polynomial /, then 
/(x) — A is SOS if and only if there exists an s x ,s matrix Q ^ for which /(x) — A = 
V(/(x)^2Vf/('^)- Now we can compare the coefficients on both sides of the latter 
equation. This leads to the SOS -program 



St. /o - A = go.o 

T.p.y.p+y=aQp.r = fa 

In a similar vein. Theorem [2T| and Proposition [3] jointly imply that we can use 
semi-definite optimization to provide a sequence of semi-definite relaxations for 



the optimization problem (21 1. Assuming that the set K satisfies the assumption of 



Theorem 21 then if /(x) - 7 > for all x e K, then 

/(x)-7 = po(x) + £p,(xk,W, (23) 



where pi(x), / = 0, 1, . . . ,m have a sum of squares representation. Theorem 21 does 
not specify the degree of the polynomials /?, (x). Thus, we select a bound 2d on the 
degree of the polynomials /?,(x), and we apply Proposition[3]to each of the polyno- 
mials Pi{x), that is, ^((x) is a sum of squares if and only if pi{x) — Vdi'x)'^ QiVd{x.) 
with Qi ^0, i = 0,1, . . . ,m. Substituting to Eq. ([23]l, we obtain that 7,2,, / ~ 



0,1,..., m, satisfy linear equations that we denote as L{y, Qo,Qi,... , Qm) = 0. Thus, 



we can find a lower bound to problem (20 1 by solving the semi-definite optimization 
problem 

= max 7 

St. L{y,QQ,Qu...,Q„^)=0, (24) 
Q,>0, i = 0,l,...,m. 



Problem ( |24] l involves semi-definite optimization over m + 1 s x s matrices. From 
the construction we get the relation Zj < Z*. It turns out that as d increases, Z^/ 
converges to Z*. Moreover, for binary optimization, there exists a finite d for which 



Problem ( 24 1 provides a systematic way to find convergent semi-definite relax- 
ations to problem (|20]l. While the approach is both general (it applies to very general 
nonconvex problems including nonlinear mixed-integer optimization problems) and 
insightful from a theoretical point of view, it is only practical for values of li = 1,2, 
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as large scale semi-definite optimization problems cannot be solved in practice. In 
many situations, however, Z\ or Zi provide strong bounds. Let us consider an exam- 
ple. 

Example 2. Let us minimize /(xi ,^2) = 2x\ + 2x\x2 — x\x\ + 5x\ over K^. We at- 
tempt to write 

f{xi 5X2) = 2x\ + 2x[X2 — XjX2 -|- 5x2 





T 


qn qi2 qu 




( -A ^ 


X2 




q 12 qii qii 






\x1x2 ) 




_qu q23 ?33_ 




^XlX2 / 



= q\ \x\ + ^22X2 + {qn + 2^l2)xiX2 + 2^i3XjX2 + 2^23x1X2- 

In order to have an identity, we obtain 

^11 =2, ^22 = 5, ^33 -h 2^12 = -1, 2^13 =2, ^23 = 0. 

Using semi-definite optimization, we find a particular solution such that 2 ^ is 
given by 

^ " 20' 



Q 



2 -3 1 
-3 5 
1 05 



HH' 



H. 



1 



-3 1 
1 3 



It follows that /(xi ,X2) = J (2xj — Sxj +xiX2)^ + j (x2 + 3x1x2)^, and thus the opti- 
mal solution value is 7* = and the optimal solution is x| = Xj = 0. 



4.3 Quadratic functions 



In this section, we focus on instances of polynomial programming where the func- 
tions are all quadratic. The specific form of the mixed-integer quadratically con- 
strained programming problem that we consider is 



min ^o(x) 
s.t. ^(x) < 

1<X<U 

Xi G M 
x,-gZ 



(MIQCP[l,u]) 



for i= I, 
for i = k- 



..,k 



where 50 : ^ R and ^ : M" ^ M'" are quadratic, 1, u G Z«, and 1 < u. We denote 
the continuous relaxation by (MIQCPr[I,u]). We emphasize that we are not gener- 
ally making any convexity/concavity assumptions on the quadratic functions qt, so 
when we do require any such assumptions we will state so exphcitly. 
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Of course one can write a binary constraint y, e {0,1} as the (nonconvex) 
quadratic inequality 3',(1 — y,) < in the bound-constrained variable < < 1. 
So, in this way, the case of binary variables y, can be seen as the special case 
of (MIQCP[l,u]) with no discrete variables (i.e., k — n). So, in a sense, the topic 
of mixed-binary quadratically constrained programming can be seen as a special 
case of (purely continuous) quadratically constrained programming. We are not 
saying that it is necessarily useful to do this from a computational viewpoint, but 
it makes it clear that the scope of even the purely continuous quadratic model in- 
cludes quadratic models having both binary and continuous variables, and in partic- 
ular mixed-{0, 1} linear programming. 

In addition to the natural mathematical interest in studying mixed-integer quadrat- 
ically constrained programming, there is a wealth of applications that have mo- 
tivated the development of practical approaches; for example: Trimloss problems 
(see |84j, for example), portfolio optimization (see 1261 . for example), Max-Cut 
and other binary quadratic models (see 1108111091 and the references therein). 

In the remainder of this section, we describe some recent work on practical com- 
putational approaches to nonconvex quadratic optimization models. Rather than at- 
tempt a detailed survey, our goal is to present a few recent and promising tech- 
niques. One could regard these techniques as belonging more to the field of global 
optimization, but in Section |5] we present material on global optimization aimed at 
more general unstructured nonlinear integer programming problems. 

4.3.1 Disjunctive programming 

It is not surprising that integer variables in a mixed-integer nonlinear program can 
be treated with disjunctive programming 17J|9J. A corresponding branch-and-cut 
method was first described in II 1231 in the context of 0/1 mixed convex program- 
ming. 

Here we describe an intriguing result from 11 14111 15l ll 16\. which shows that one 
can also make useful disjunctions from nonconvex quadratic functions in a mixed- 
integer quadratic ally-constrained programming problem or even in a purely con- 
tinuous quadratic ally-constrained programming problem. The starting point for this 
approach is that we can take a quadratic form x^A,x in x e M", and rewrite it via an 
extended formulation as the linear form (A;,X), using the matrix variable X e M"^", 
and the nonlinear equation X xx^. The standard approach is to relax X = xx^ to 
the convex inequality X >z xx^. But the approach of 11 14111 15111161 involves work- 
ing with the nonconvex inequality X ^ xx^. This basic idea is as follows. Let v G M" 
be arbitrary (for now). We have the equation 

(vvT,X) = (vv^.xx"^) = (v'^x)2 , 

which we relax as the concave inequality 



(vTx)2>(vvT,X). 
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If we have a point (x,X) that satisfies the convex inequality X y xx^ , but for which 
X 7^ xx^ , then it is the case that X — xx^ has a positive eigenvalue X . Let v denote 
a unit-length eigenvector belonging to X . Then 



So, A > if and only if (v^x)^ < (vv^,X) . That is, every positive eigenvalue 
of X — xx^ yields an inequality of the form (Q) that is violated by (x,X) . Next, 
we make a disjunction on this violated nonconvex inequality (£2). First, we choose 
a suitable polyhedral relaxation ^ of the feasible region, and we let [t7l,T7j/] be 
the range of v^x as (x,X) varies over the relaxation Next, we choose a value 
G (tJl,T]£/) (e.g., the midpoint), and we get the polyhedral disjunction; 



Notice that the second part of the first (resp., second) half of the disjunction corre- 
sponds to a secant inequality over the interval between the point 6 and the lower 
(resp., upper) bound for v^x. Finally, we use the linear-programming technology 
of ordinary disjunctive programming to separate, via a linear inequality, the point 
(x,X) from the convex closure of the two halves of the disjunction. Details and 
extensions of this idea appear in II114II115TI116L 

4.3.2 Branch and cut 

A branch-and-cut scheme for optimization of a nonconvex quadratic form over a box 
was recently developed by Vandenbussche and Nemhauser 11331 11321 . They use 
a formulation of Balas via linear programming with complementarity conditions, 
based on the necessary optimality conditions of continuous quadratic programming 
(see (HI). Specifically, they consider the problem 



A=A||v||^ 

= (vv"^,X-xx"r). 




I7l(v) < v^x < 

• (vTx)(77i(v) + 0) 




or 




min jX^Qx + c^x 

s.t. xe [0,1]", 



(BoxQP[e,c]) 



where 2 is an n x n symmetric, non positive semi-definite matrix, and c £ M". The 
KKT necessary optimality conditions for (BoxQP[2,c]) are 
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y-ex-z = c, (25) 

yT(l-x) = (26) 

z^x = (27) 

X e [0, 1]" (28) 

y, z e M'l . (29) 

Vandenbussche and Nemhauser, appealing to a result of Balas, define ^(2,c) as 
the polyhedron defined as the convex hull of solutions to ( 25p9 1, and they work 
with the reformulation of (BoxQP) as the linear program 

min ic''"x+ il''"y , , 

s.t. (x,y,z)l^(e,c). (^^i^^[e,c]) 

The main tactic of Vandenbussche and Nemhauser is to develop cutting planes for 

^(e,c). 

Burer and Vandenbussche pursue a similar direction, but they allow general poly- 
hedral constraints and employ semi -definite relaxations 



4.3.3 Branch and bound 

Linderoth also looks at quadratically-constrained programs that are not convex ||9T1l . 
He develops a novel method for repeatedly partitioning the continuous feasible re- 
gion into the Cartesian product of triangles and rectangles. What is particularly inter- 
esting is that to do this effectively, Linderoth develops convex envelopes of bilinear 
functions over rectangles and triangles (also see Anstreicher and Burer's paper |6 1), 
and then he demonstrates that these envelopes involve hyperbolic constraints which 
can be reformulated as the second-order cone constraints. It is interesting to com- 
pare this with the similar use of second-order cone constraints for convex quadratics 
(see Section [3.3.2| i. 

One can view the technique of Linderoth as being a specialized "Spatial Branch- 
and-Bound Algorithm." In Section|5]we will describe the Spatial Branch-and-Bound 
Algorithm for global optimization in its full generality. 



5 Global optimization 

In the present section we take up the subject of global optimization of rather gen- 
eral nonlinear functions. This is an enormous subject, and so we will point to just a 



couple of directions that we view as promising. On the practical side, in Section 5.1 



we describe the Spatial Branch-and-Bound Algorithm which is one of the most suc- 



cessful computational approaches in this area. In Section 5.2 from the viewpoint of 



complexity theory, with a goal of trying to elucidate the boundary between tractable 
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and intractable, we describe some very recent work on global optimization of a very 
general class of nonlinear functions over an independence system. 



5.1 Spatial Branch-and-Bound 

In this section we address methods for global optimization of rather general mixed- 
integer nonlinear programs having non-convex relaxations. Again, to have any hope 
at all, we assume that the variables are bounded. There is a very large body of 
work on solution techniques in this space. We will not attempt to make any kind of 
detailed survey. Rather we will describe one very successful methodology, namely 
the Spatial Branch-and-Bound Algorithm. We refer to Ill22lll27 ll and the references 
therein. 

The Spatial Branch-and-Bound Algorithm for mixed-integer nonlinear program- 
ming has many similarities to the ordinary branch-and-bound employed for the so- 
lution of mixed-integer linear programs, but there are many additional wrinkles. 
Moreover, the techniques can be integrated. In what follows, we will concentrate on 
how continuous nonlinearities are handled. We leave it to the reader to see how these 
techniques would be integrated with the associated techniques for mixed-integer lin- 
ear programs. 

One main difference with the mixed-integer Unear case is that all nonlinear func- 
tions in a problem instance are symbolically and recursively decomposed via simple 
operators, until we arrive at simple functions. The simple operators should be in a 
limited library. For example: sum, product, quotient, exponentiation, power, loga- 
rithm, sine, cosine, absolute value. Such a decomposition is usually represented via 
a collection of rooted directed acyclic graphs. At each root is a nonlinear function 
occurring in the problem formulation. Leaves are constants, affine functions and 
atomic variables. Each non-leaf node is thought of as an auxiliary variable and also 
as representing a simple operator, and its children are the arguments of that operator. 

An inequality constraint in the problem formulation can be thought of as a bound- 
ing interval on a root. In addition, the objective function is associated with a root, 
and so lower and upper bounds on the optimal objective value can also be thought 
of as a bounding interval on a root. Simple bounds on a variable in the problem 
formulation can be thought of as a bounding interval on a leaf. In this way, we have 
an extended-variable reformulation of the given problem. 

Bounds are propagated up and down each such rooted directed acycUc graph via 
interval arithmetic and a library of convex envelopes or at least linear convex relax- 
ations of the graphs of simple nonlinear operators acting on one or two variables on 
simple domains (intervals for univariate operators and simple polygons for bivariate 
operators). So, in this way, we have a now tractable convex or even linear relaxation 
of the extended-variable reformulation, and this is used to get a lower bound on the 
minimum objective value. 

The deficiency in our relaxation is localized to the graphs of simple functions that 
are only approximated by convex sets. We can seek to improve bounds by branching 
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on the interval for a variable and reconvexifying on the subproblems. For example, 
we may have a variable v that is a convex function / in a variable w on an inter- 
val [l,u]. Then the convex envelope of the graph G[/;m] :— {{v,w) : v = f{w)} is 
precisely 

G[/,«]:=|(v,w) : /(w) < v </(/)+ (^-^^^^^^^ (w- Z)} • 

We may find that at a solution of the relaxation, the values of the variables (m,v), 
say (v,M>) € G[Z,m], are far from G[l,u]. In such a case, we can branch by choosing 
a point b E [l,u] (perhaps at or near v), and forming two subproblems in which 
the bounding interval for v is amended as [l,b] in one and [b,u] in the other The 
value in branching on such a continuous variable is that we now reconvexify on 
the subintervals, effectively replacing G[l, u] with the smaller set G[l,b](J G[b,u]. In 
particular, if we did choose b~v, then (v,m>) ^ G[l ,b](JG[l,b], and so the algorithm 
makes some progress. We note that a lot of work has gone into good branching 
strategies (see [14| for example). 

Finally, a good Spatial Branch-and-Bound procedure should have an effective 
strategy for finding good feasible solutions, so as to improve the objective upper 
bound (for minimization problems). A good rudimentary strategy is to take the so- 
lution of a relaxation as a starting point for a continuous nonlinear-programming 
solver aimed at finding a locally-optimal solution of the continuous relaxation (of 
either the original or extended-variable formulation). Then if a feasible solution to 
this relaxation is obtained and if it happens to have integer values for the appro- 
priate variables, then we have an opportunity to update the objective value upper 
bound. Alternatively, one can use a solver aimed mainly at mixed-integer nonlinear 
programs having convex relaxation as a heuristic also from such a starting point. 
In fact, the Branch-and-Bound Algorithm in Bonmin has options aimed at giving 
good solutions from such a starting point, even for non-convex problems. 

The Spatial Branch-and-Bound Algorithm relies on the rapid and tight convex- 
ification of simple functions on simple domains. Therefore, considerable work has 
gone into developing closed-form expressions for such envelopes. This type of work 
has paralleled some research in mixed-integer linear programming that has focused 
on determining convex hulls for simple constraints. Useful results include: univari- 
ate functions ll2l 11211 135]| . univariate monomials of odd degree 1891 l90ll , bilinear 
functions ||5]|97l, trilinear functions ||98ll , so-called (n— l)-convex functions ||73]| . 
and fractional terms B126L Further relevant work includes algorithms exploiting 
variable transformations and appropriate convex envelopes and relaxations. For ex- 
ample, for the case of "signomials" (i.e., terms of the form oqx"^ ■ ■ ■ x^]" , with 
at G M), see lITOl . 

We do not make any attempt to exhaustively review available software for global 
optimization. Rather we just mention that state-of-the-art codes implementing a 
Spatial Branch-and-Bound Algorithm include Baron Ml 12II1271I128I and the new 
open-source code Couenne fT4l|. 
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5.2 Boundary cases of complexity 

Now, we shift our attention back to the viewpoint of complexity theory. Our goal is 
to sample a bit of the recent work that is aimed at revealing the boundary between 
tractable and intractable instances of nonlinear discrete optimization problems. We 
describe some very recent work on global optimization of a very general class of 
nonlinear functions over an independence system (see |85|). Other work in this vein 
includes 061 [TTl . 

Specifically, we consider the problem of optimizing a nonlinear objective func- 
tion over a weighted independence system presented by a linear-optimization oracle. 
While this problem is generally intractable, we are able to provide a polynomial- 
time algorithm that determines an "r-best" solution for nonlinear functions of the 
total weight of an independent set, where r is a constant that depends on certain 
Frobenius numbers of the individual weights and is independent of the size of the 
ground set. 

An independence system is a nonempty set of vectors S C {0, 1}" with the prop- 
erty that X g {0, 1}", X < y e 5 implies x e 5. The general nonlinear optimization 
problem over a multiply-weighted independence system is as follows. Given an in- 
dependence system S C {0,1}", weight vectors wi,...,Wf/ G Z", and a function 
/ : Z'' ^ M, find x e 5 minimizing the objective /(w^x, . . . ,wjx). 

The representation of the objective in the above composite form has several ad- 
vantages. First, for c/ > 1, it can naturally be interpreted as multi-criteria optimiza- 
tion: the d given weight vectors wi , . . . ,W(/ represent d different criteria, where the 
value of X G 5 under criterion / is its /-th total weight w^x and the objective is 
to minimize the "balancing" /(w^^x, . . . , wjx) of the d given criteria by the given 
function /. Second, it allows us to classify nonlinear optimization problems into a 
hierarchy of increasing generality and complexity: at the bottom lies standard linear 
optimization, recovered with d = \ and / the identity on Z; and at the top lies the 
problem of maximizing an arbitrary function, which is typically intractable, arising 
with d = n and w, = 1, the /-th standard unit vector in Z" for all /. 

The computational complexity of the problem depends on the number d of 
weight vectors, on the weights w, j, on the type of function / and its presentation, 
and on the type of independence system S and its presentation. For example, when 
5 is a matroid, the problem can be solved in polynomial time for any fixed d, any 
{0, 1, . . . ,p}-valued weights Wij with p fixed, and any function / presented by a 
comparison oracle, even when S is presented by a mere membership oracle, see 
iri6J . Also, for example, when S consists of the matchings in a given bipartite graph 
G, the problem can be solved in polynomial time for any fixed d, any weights Wij 
presented in unary, and any concave function /, see flSl; but on the other hand, for 
convex f, already with fixed d — 2 and {0, 1} -valued weights Wij, the problem in- 
cludes as a special case the exact matching problem whose complexity is long open 
11100111021 . 

In view of the difficulty of the problem already for d ~2,we take a first step and 
concentrate on nonlinear optimization over a (singly) weighted independence sys- 
tem, that is, with d = I, single weight vector w = (wi , . . . , w„), and univariate func- 
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tion / : Z ^ M. The function / can be arbitrary and is presented by a comparison 
oracle that, queried on x,y e Z, asserts whether or not /(x) < /(y). The weights wj 
take on values in a /9-tuple a — {ai,... ,ap) of positive integers. Without loss of gen- 
erality we assume that a — {ai, . . . , ap) is primitive, by which we mean that the o; are 
distinct positive integers whose greatest common divisor gcd(a) :— gcd(fli , . . . ,flp) 
is 1. The independence system S is presented by a linear-optimization oracle that, 
queried on vector c G Z", returns an element x e 5 that maximizes the linear function 
c^x = Y!]=i '^i^j- It turns out that this problem is akeady quite intriguing, and so we 
settle for an approximative solution in the following sense, that is interesting in its 
own right. For a nonnegative integer r, we say that x* e 5 is an r-best solution to the 
optimization problem over S if there are at most r better objective values attained by 
feasible solutions. In particular, a 0-best solution is optimal. Recall that the Frobe- 
nius number of a primitive a is the largest integer F(a) that is not expressible as a 
nonnegative integer combination of the a,-. We prove the following theorem. 

Theorem 22. For every primitive p-tuple a — (ai, . . . ,flp), there is a constant r{a) 
and an algorithm that, given any independence system S C {0, 1}" presented by a 
linear-optimization oracle, weight vector w G {ai , . . . , ap}", and function / : Z ^ M 
presented by a comparison oracle, provides an r(a)-best solution to the nonlinear 
problem min{/(w^x) : x G 5}, in time polynomial in n. Moreover: 

1. If ai divides fl,+ i for / = 1, . . . ,p — 1, then the algorithm provides an optimal 
solution. 

2. For p = 2, that is, for a~ (oi , a2)> the algorithm provides an F{a)-best solution. 

In fact, we give an explicit upper bound on r(a) in terms of the Frobenius num- 
bers of certain subtuples derived from a. An interesting special case is that of 
a = (2, 3). Because F(2, 3) = 1, the solution provided by our algorithm in that case 
is either optimal or second best. 



The proof of Theorem 22 is pretty technical, so we only outline the main ideas. 
Below we present a naive solution strategy that does not directly lead to a good ap- 
proximation. However, this naive approach is used as a basic building block. One 
partitions the independence system into suitable pieces, to each of which a suitable 
refinement of the naive strategy is applied separately. Considering the monoid gen- 
erated by {fli , . . . ,flp} allows one to show that the refined naive strategy applied to 
each piece gives a good approximation within that piece. In this way, the approx- 
imation quality r(a) can be bounded as follows, establishing a proof to Theorem 

EH 

Lemma 6. Let a = (oi , . . . , ap) be any primitive p-tuple. Then the following hold: 

1. An upper bound on ria) is given by r(a) < (2max(a))''. 

2. For divisible a, we have r(a) = 0. 

3. For p — 2, that is, for a = (oi ,02), we have r(a) = F(a). 

Before we continue, let us fix some notation. The indicator of a subset J CN is 
the vector Ij := Ljey 1; € {Oi 1 }", so that supp(l j) = J. Unless otherwise specified. 
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X denotes an element of {0,1}" and A,T,V denote elements of Z'^. Throughout, 
a = (ai,...,flp) is a primitive /?-tuple. We will be working with weights taking 
values in a, that is, vectors w G {oi , . . . ,flp}". With such a weight vector w being 
clear from the context, we let A^,- := {j £ N : Wj = a,} for / = l,...,p, so that 
N = l+lf^i A^;. For X e {0, 1}" we let A,(x) |supp(x) n A^,| for / = 1, . . . ,p, and 
A(x) := (Ai (x), . . . , Ap(x)), so that w'''x = A(x)^a. For integers z,s £Z and a set of 
integers Z C Z, we define z + sZ := {z + sx : x G Z}. 

Let us now present the naive strategy to solve the univariate nonlinear problem 
min{/(w^x) : x G S}. Consider a set 5 C {0, 1}", weight vector w £ {fli, . . . ,flp}", 
and function / : Z ^ M presented by a comparison oracle. Define the image of S 
under w to be the set of values w^x taken by elements of S; we denote it by w • S. 

We point out the following simple observation. 

Proposition 4. A necessary condition for any algorithm to find an r-best solution to 
the problem min{/(w^x) : x e 5}, where the function f is presented by a compar- 
ison oracle only, is that it computes all but at most r values of the image v/ ■ S of S 
under w. 

Note that this necessary condition is also sufficient for computing the objective value 
/(w^x*) of an r-best solution, but not for computing an actual r-best solution x* e 
S, which may be harder Any point x attaining max{w^x ; x e 5} provides an 
approximation of the image given by 

{w^x : x<x}Cw5C{0,l,...,w"^x} . (30) 

This suggests the following natural naive strategy for finding an approximative so- 
lution to the optimization problem over an independence system S that is presented 
by a linear-optimization oracle. 

Algoritlim 4 (Naive Strategy) 

input Independence system S C {0, 1}" presented by a linear-optimization oracle, 

/ ; Z — > M presented by a comparison oracle, and w S {ai , . . . ,ap}" 

obtain x attaining maxjw^x : x G 5} using the linear- optimization oracle for S 

output X* as one attaining min{/(w^x) : x < x} using the algorithm of Lemma^ 

below. 

Unfortunately, as the next example shows, the number of values of the image that 



are missing from the approximating set on the left-hand side of equation ( 30 1 cannot 
generally be bounded by any constant. So by Proposition |4j this strategy cannot be 
used as is to obtain a provably good approximation. 

Example 3. Let a (1,2), n 4m, y ^4=1 1/' ^ := Lf=2m+i 1'' ^ := y + 2z, 
that is, 

y = (1,...,1,0,...,0), z = (0,...,0,1,...,1), w = (1,...,1,2,...,2), 



define / on Z by 
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^odd; 
[ 2m, k even, 

and let S be the independence system 

5 := {x e {0, 1}" : X < y} U {x e {0, 1}" : X < z}. 

Then the unique optimal solution of the linear-objective problem maxjw^x ; x e 5} 
is X : = z, with w^x = Am, and therefore 

{w^x : X < x} = {2i : i — 0,1, ... , 2m} and 

w-S' = {/ : i = 0,1, ... ,2m} U {2i : i — 0,1, ... ,2m}. 

So all m odd values (i.e., 1,3, .. . ,2m — 1) in the image w • S are missing from the 
approximating set {w^x : x < x} on the left-hand side of ( 30 1, and x* attaining 
min{/(w^x) : x < x} output by the above strategy has objective value /(w^x*) = 
2m, while there are m — ^ better objective values (i.e., 1,3, .. . ,2m — 1) attainable 
by feasible points (e.g., Y!i=i I/' for k= 1,3, .. . ,2m — 1). 

Nonetheless, a more sophisticated refinement of the naive strategy, applied re- 
peatedly to several suitably chosen subsets of S rather than S itself, will lead to a 
good approximation. Note that the naive strategy can be efficiently implemented as 
follows. 

Lemma 7. For every fixed p-tuple a, there is a polynomial-time algorithm that, 
given univariate function / : Z — > M presented by a comparison oracle, weight vec- 
for w €{«!,... ,flp}", an^ix G {0, 1}", io/vei min{/(w^x) : x < x}. 

Proof. Consider the following algorithm: 

Algorithm 5 input function / : Z ^ M presented by a comparison oracle, w S 

{oi, . . .,ap}" andx G {0, 1}" 

let Af, := {j : Wj = a,} and T; A,(x) = |supp(x) nM|, / — l,...,p 

For every choice of V = (Vi , . . . , Vp) < (ti , . . . , Tp) = T 

determine some Xy < x with A,(xv) = |supp(xv) nA^,| — Vi, i — I, . . . ,p 

output X* as one minimizing /(w^x) among the Xy by using the comparison 

oracle of f. 

As the value w^x depends only on the cardinalities |supp(x) nA^,|, i— 1,. .. ,p, 
it is clear that 

{w^x : X < x} = {w^Xv : V < t}. 

Clearly, for each choice V < T it is easy to determine some Xy < x by zeroing out 
suitable entries of x. The number of choices V < T and hence of loop iterations and 
comparison-oracle queries of / to determine x* is 



n(T,+i) < («+!)". 
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Table 1 Computational complexity and algorithms for nonlinear integer optimization. 



Objective 
function 



Linear 



Constraints 
Convex Polynomial 



Arbitrary Polynomial 



Linear Polynomial-time in fixed 

dimension: 

- Lenstra's algorithm f87 1 

- Generalized basis re- 
duction, Lovasz-Scarf 
|94| 

- Short rational generating 
functions, Barvinok 1121 

Convex max Polynomial-time in fixed 
(§|2) dimension: Cook et al. 138) 



Polynomial-time in fixed 
dimension: 

Lenstra-type algorithms 
(§|3l) 



Khachiyan-Porkolab 
Heinz 1661 



Incomputable: 

Hubert's 10th problem, 
Matiyasevich |95| (§[l), 
even for: 

- quadratic constraints, 
Jeroslow |75 1 

- fixed dimension 10, 
Matiyasevich 1761 



Incomputable (§ 1 



Convex min Polynomial-time in fixed dimension: Lenstra-type al- Incomputable (j 
(§|3j gorithms: Khachiyan-Porkolab 1 79 1, Heinz |66] (§ [3ll) 

Arbitrary NP-hard, inapproximable, even for quadratic forms Incomputable ([ 

Polynomial over hypercubes: MAX-CUT, Hastad |65 | (§[!} 

(§1^ NP-hard, even for fixed dimension 2, degree 4 

FPTAS in fixed dimen- 
sion: Short rational gener- 
ating functions, De Loera 
et al. tlJJ (§[4l) 



6 Conclusions 



In this chapter, we hope to have succeeded in reviewing mixed-integer nonlinear 
programming from two important viewpoints. 

We have reviewed the computational complexity of several important classes 
of mixed-integer nonlinear programs. Some of the negative complexity results (in- 
computability, NP-hardness) that appeared in Section[T]have been supplemented by 
polynomiality or approximability results in fixed dimension. Table [T] gives a sum- 
mary. In addition to that, and not shown in the table, we have explored the bound- 
ary between tractable and intractable problems, by highlighting interesting cases in 
varying dimension where still polynomiality results can be obtained. 

Additionally, we have reviewed a selection of practical algorithms that seem to 
have the greatest potential from today's point of view. Many of these algorithms, 
at their core, are aimed at integer convex minimization. Here we have nonlinear 
branch-and-bound, outer approximation, the Quesada-Grossman algorithm, hybrid 
algorithms, and generalized Benders decomposition. As we have reported, such ap- 



52 



Raymond Hemmecke, Matthias Koppe, Jon Lee and Robert Weismantel 



proaches can be specialized and enhanced for problems with SDP constraints, SOCP 
constraints, and (convex) quadratics. For integer polynomial programming (without 
convexity assumptions), the toolbox of Positivstellensatze and SOS programming is 
available. For the case of quadratics (without convexity assumptions), specialized 
versions of disjunctive programming, branch-and-cut, and branch-and-bound have 
been devised. Finally, for general global optimization, spatial branch-and-bound is 
available as a technique, which relies heavily on convexification methods. 

It is our hope that, by presenting these two viewpoints to the interested reader, 
this chapter will help to create a synergy between both viewpoints in the near future. 
May this lead to a better understanding of the field, and to much better algorithms 
than what we have today! 
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